simsopt.mhd package
- class simsopt.mhd.Boozer(equil: Vmec, mpol: int = 32, ntor: int = 32)
Bases:
Optimizable
This class handles the transformation to Boozer coordinates.
A Boozer instance maintains a set “s”, which is a registry of the surfaces on which other objects want Boozer-coordinate data. When the run() method is called, the Boozer transformation is carried out on all these surfaces. The registry can be cleared at any time by setting the s attribute to {}.
- 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.
- register(s: Union[float, Iterable[float]]) None
This function is called by objects that depend on this Boozer object, to indicate that they will want Boozer data on the given set of surfaces.
- Parameters
s – 1 or more surfaces on which Boozer data will be requested.
- run()
Run booz_xform on all the surfaces that have been registered.
- class simsopt.mhd.IotaTargetMetric(vmec, iota_target, adjoint_epsilon=0.1)
Bases:
Optimizable
IotaTargetMetric is a class that computes a metric quantifying the deviation of the rotational transform \(\iota\) in from a prescribed target profile in a Vmec equilibrium:
\[J = \frac{1}{2} \int ds \, (\iota - \iota_{target})^2\]where the integral is over the normalized toroidal flux \(s\), and the function \(\iota_{target}(s)\) corresponds to the argument
iota_target
. This class also can compute the derivatives of \(J\) using an adjoint method.- Parameters
vmec – instance of Vmec
iota_target – function handle which takes a single argument, s, the normalized toroidal flux, and returns the target rotational transform.
adjoint_epsilon – sets the amplitude of the toroidal current perturbation required for the adjoint solve.
- J()
Computes the quantity \(J\) described in the class definition.
- dJ()
Computes derivatives of \(J\) with respect to surface parameters using an adjoint method.
- shape_gradient()
Computes the shape gradient of the quantity \(J\) described in the class definition. For a perturbation to the surface \(\delta \vec{x}\), the resulting perturbation to the objective function is
\[\delta J(\delta \vec{x}) = \int d^2 x \, G \delta \vec{x} \cdot \vec{n}\]where the integral is over the VMEC boundary surface, \(G\) is the shape gradient, and \(\vec{n}\) is the unit normal.
- Returns
2d array of size (numquadpoints_phi,numquadpoints_theta)
- Return type
\(G\)
- class simsopt.mhd.IotaWeighted(vmec, weight_function, adjoint_epsilon=0.1)
Bases:
Optimizable
Computes a weighted average of the rotational transform for a VMEC configuration. The quantity computed is defined by
\[J = \frac{ \int ds \, \iota(s) w(s)} { \int ds \, w(s)}\]where \(w(s)\) is a prescribed weight function, corresponding to the argument
weight_function
. This class also can compute the derivatives of \(J\) using an adjoint method.- Parameters
vmec – instance of Vmec
weight_function – function handle which takes a single argument, s, the normalized toroidal flux
adjoint_epsilon – sets the amplitude of the toroidal current perturbation required for the adjoint solve.
- J()
Computes the quantity \(J\) described in the class definition.
- dJ()
Computes derivatives of \(J\) with respect to surface parameters using an adjoint method.
- shape_gradient()
Computes the shape gradient of the quantity \(J\) described in the class definition. For a perturbation to the surface \(\delta \vec{x}\), the resulting perturbation to the objective function is
\[\delta J(\delta \vec{x}) = \int d^2 x \, G \delta \vec{x} \cdot \vec{n}\]where the integral is over the VMEC boundary surface, \(G\) is the shape gradient, and \(\vec{n}\) is the unit normal.
- Returns
2d array of size (numquadpoints_phi,numquadpoints_theta)
- Return type
\(G\)
- class simsopt.mhd.Profile(*args, **kwargs)
Bases:
Optimizable
Base class for radial profiles. This class should not be used directly - use subclasses instead.
- __call__(*args, **kwargs)
Shortcut for calling f(s)
- plot(ax=None, show=True, n=100)
Plot the profile using matplotlib.
- Parameters
ax – The axis object on which to plot. If
None
, a new figure will be created.show – Whether to call matplotlib’s
show()
function.n – The number of grid points in s to show.
- class simsopt.mhd.ProfilePolynomial(data)
Bases:
Profile
A profile described by a polynomial in the normalized toroidal flux s. The polynomial coefficients are dofs that are fixed by default.
- Parameters
data – 1D numpy array of the polynomial coefficients. The first coefficient is the constant term, the next coefficient is the linear term, etc.
- dfds(s)
Return the d/ds derivative of the profile at specified points in s.
- f(s)
Return the value of the profile at specified points in s.
- class simsopt.mhd.ProfilePressure(*args)
Bases:
Profile
A Profile \(f(s)\) which is determined by other profiles \(f_j(s)\) as follows:
\[f(s) = \sum_j f_{2j}(s) f_{2j+1}(s).\]This is useful for creating a pressure profile in terms of density and temperature profiles, with any number of species. Typical usage is as follows:
ne = ProfilePolynomial(1.0e20 * np.array([1.0, 0.0, 0.0, 0.0, -1.0])) Te = ProfilePolynomial(8.0e3 * np.array([1.0, -1.0])) nH = ne TH = ProfilePolynomial(7.0e3 * np.array([1.0, -1.0])) pressure = ProfilePressure(ne, Te, nH, TH)
This class does not have any optimizable dofs.
- Parameters
args – An even number of Profile objects.
- dfds(s)
Return the d/ds derivative of the profile at specified points in s.
- f(s)
Return the value of the profile at specified points in s.
- class simsopt.mhd.ProfileScaled(base, scalefac)
Bases:
Profile
A Profile which is equivalent to another Profile object but scaled by a constant. This constant is an optimizable dof, which is fixed by default.
- Parameters
base – A Profile object to scale
scalefac – A number by which the base profile will be scaled.
- dfds(s)
Return the d/ds derivative of the profile at specified points in s.
- f(s)
Return the value of the profile at specified points in s.
- class simsopt.mhd.ProfileSpline(s, f, degree=3)
Bases:
Profile
A Profile that uses spline interpolation via scipy.interpolate.InterpolatedUnivariateSpline
The
f
data are optimizable dofs, which are fixed by default.- Parameters
s – A 1d array with the x coordinates for the spline.
f – A 1d array with the y coordinates for the spline.
degree – The polynomial degree of the spline. Must be in
[1, 2, 3, 4, 5]
.
- dfds(s)
Return the d/ds derivative of the profile at specified points in s.
- f(s)
Return the value of the profile at specified points in s.
- resample(new_s, degree=None)
Return a new
ProfileSpline
object that has different grid points (spline nodes). The data from the old s grid will be interpolated onto the new s grid.- Parameters
new_s – A 1d array representing the x coordinates of the new
ProfileSpline
.degree – The polynomial degree used for the new
ProfileSpline
object. IfNone
, the degree of the originalProfileSpline
will be used.
- Returns
A new
ProfileSpline
object, in which the data have been resampled ontonew_s
.
- class simsopt.mhd.Quasisymmetry(boozer: Boozer, s: Union[float, Iterable[float]], helicity_m: int, helicity_n: int, normalization: str = 'B00', weight: str = 'even')
Bases:
Optimizable
This class is used to compute the departure from quasisymmetry on a given flux surface based on the Boozer spectrum.
- Parameters
boozer – A Boozer object on which the calculation will be based.
s – The normalized toroidal magnetic flux for the flux surface to analyze. Should be in the range [0, 1].
helicity_m – The poloidal mode number of the symmetry you want to achive. The departure from symmetry
B(helicity_m * theta - nfp * helicity_n * zeta)
will be reported.helicity_n – The toroidal mode number of the symmetry you want to achieve. The departure from symmetry
B(helicity_m * theta - nfp * helicity_n * zeta)
will be reported.normalization – A uniform normalization applied to all bmnc harmonics. If
"B00"
, the symmetry-breaking modes will be divided by the m=n=0 mode amplitude on the same surface. If"symmetric"
, the symmetry-breaking modes will be divided by the square root of the sum of the squares of all the symmetric modes on the same surface. This is the normalization used in stellopt.weight – An option for a m- or n-dependent weight to be applied to the bmnc amplitudes.
- J() Union[Sequence[Real], NDArray[Any, float64]]
Carry out the calculation of the quasisymmetry error.
- Returns
1D numpy array listing all the normalized mode amplitudes of symmetry-breaking Fourier modes of
|B|
.
- __init__(boozer: Boozer, s: Union[float, Iterable[float]], helicity_m: int, helicity_n: int, normalization: str = 'B00', weight: str = 'even') None
Constructor
- recompute_bell(parent=None)
Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.
This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.
Need to be implemented by classes that provide a dof_setter for external handling of DOFs.
- return_fn_map: Dict[str, Callable] = {'J': <function Quasisymmetry.J>}
- class simsopt.mhd.QuasisymmetryRatioResidual(vmec: Vmec, surfaces: Union[float, Sequence[Real], NDArray[Any, float64]], helicity_m: int = 1, helicity_n: int = 0, weights: Optional[Union[Sequence[Real], NDArray[Any, float64]]] = None, ntheta: int = 63, nphi: int = 64)
Bases:
Optimizable
This class provides a measure of the deviation from quasisymmetry, one that can be computed without Boozer coordinates. This metric is based on the fact that for quasisymmetry, the ratio
\[(\vec{B}\times\nabla B \cdot\nabla\psi) / (\vec{B} \cdot\nabla B)\]is constant on flux surfaces.
Specifically, this class represents the objective function
\[f = \sum_{s_j} w_j \left< \left[ \frac{1}{B^3} \left( (N - \iota M)\vec{B}\times\nabla B\cdot\nabla\psi - (MG+NI)\vec{B}\cdot\nabla B \right) \right]^2 \right>\]where the sum is over a set of flux surfaces with normalized toroidal flux \(s_j\), the coefficients \(w_j\) are user-supplied weights, \(\left< \ldots \right>\) denotes a flux surface average, \(G(s)\) is \(\mu_0/(2\pi)\) times the poloidal current outside the surface, \(I(s)\) is \(\mu_0/(2\pi)\) times the toroidal current inside the surface, \(\mu_0\) is the permeability of free space, \(2\pi\psi\) is the toroidal flux, and \((M,N)\) are user-supplied integers that specify the desired helicity of symmetry. If the magnetic field is quasisymmetric, so \(B=B(\psi,\chi)\) where \(\chi=M\vartheta - N\varphi\) where \((\vartheta,\varphi)\) are the poloidal and toroidal Boozer angles, then \(\vec{B}\times\nabla B\cdot\nabla\psi \to -(MG+NI)(\vec{B}\cdot\nabla\varphi)\partial B/\partial \chi\) and \(\vec{B}\cdot\nabla B \to (-N+\iota M)(\vec{B}\cdot\nabla\varphi)\partial B/\partial \chi\), implying the metric \(f\) vanishes. The flux surface average is discretized using a uniform grid in the VMEC poloidal and toroidal angles \((\theta,\phi)\). In this case \(f\) can be written as a finite sum of squares:
\[f = \sum_{s_j, \theta_j, \phi_j} R(s_j, \theta_k, \phi_{\ell})^2\]where the \(\phi_{\ell}\) grid covers a single field period. Here, each residual term is
\[R(\theta, \phi) = \sqrt{w_j \frac{n_{fp} \Delta_{\theta} \Delta_{\phi}}{V'}|\sqrt{g}|} \frac{1}{B^3} \left( (N-\iota M)\vec{B}\times\nabla B\cdot\nabla\psi - (MG+NI)\vec{B}\cdot\nabla B \right).\]Here, \(n_{fp}\) is the number of field periods, \(\Delta_{\theta}\) and \(\Delta_{\phi}\) are the spacing of grid points in the poloidal and toroidal angles, \(\sqrt{g} = 1/(\nabla s\cdot\nabla\theta \times \nabla\phi)\) is the Jacobian of the \((s,\theta,\phi)\) coordinates, and \(V' = \int_0^{2\pi} d\theta \int_0^{2\pi}d\phi |\sqrt{g}| = dV/d\psi\) where \(V\) is the volume enclosed by a flux surface.
- Parameters
vmec – A
simsopt.mhd.vmec.Vmec
object from which the quasisymmetry error will be calculated.surfaces – Value of normalized toroidal flux at which you want the quasisymmetry error evaluated, or a list of values. Each value must be in the interval [0, 1], with 0 corresponding to the magnetic axis and 1 to the VMEC plasma boundary. This parameter corresponds to \(s_j\) above.
helicity_m – Desired poloidal mode number \(M\) in the magnetic field strength \(B\), so \(B = B(s, M \vartheta - n_{fp} \hat{N} \varphi)\) where \(\vartheta\) and \(\varphi\) are Boozer angles.
helicity_n – Desired toroidal mode number \(\hat{N} = N / n_{fp}\) in the magnetic field strength \(B\), so \(B = B(s, M \vartheta - n_{fp} \hat{N} \varphi)\) where \(\vartheta\) and \(\varphi\) are Boozer angles. Note that the supplied value of
helicity_n
will be multiplied by the number of field periods \(n_{fp}\), so typicallyhelicity_n
should be +1 or -1 for quasi-helical symmetry.weights – The list of weights \(w_j\) for each flux surface. If
None
, a weight of \(w_j=1\) will be used for all surfaces.ntheta – Number of grid points in \(\theta\) used to discretize the flux surface average.
nphi – Number of grid points per field period in \(\phi\) used to discretize the flux surface average.
- compute()
Compute the quasisymmetry metric. This function returns an object that contains (as attributes) all the intermediate quantities for the calculation. Users do not need to call this function for optimization; instead the
residuals()
function can be used. However, this function can be useful if users wish to inspect the quantities going into the calculation.
- profile()
Return the quasisymmetry metric in terms of a 1D radial profile. The residuals \(R\) are squared and summed over theta and phi, but not over s. The total quasisymmetry error \(f\) returned by the
total()
function is the sum of the values in the profile returned by this function.
- residuals()
Evaluate the quasisymmetry metric in terms of a 1D numpy vector of residuals, corresponding to \(R\) in the documentation for this class. This is the function to use when forming a least-squares objective function.
- total()
Evaluate the quasisymmetry metric in terms of the scalar total \(f\).
- class simsopt.mhd.RedlGeomBoozer(booz, surfaces, helicity_n, ntheta=64, plot=False)
Bases:
Optimizable
Evaluate geometry data needed to evaluate the Redl bootstrap current formula, such as the effective fraction of trapped particles. In the approach here, Boozer coordinates are computed, and all the symmetry-breaking Bmn harmonics are discarded to obtain an effectively perfectly quasisymmetric configuration.
- Parameters
booz – An instance of
simsopt.mhd.boozer.Boozer
surfaces – A 1D array with the values of normalized toroidal flux to use for the bootstrap current calculation.
helicity_n – 0 for quasi-axisymmetry, or +/- 1 for quasi-helical symmetry. This quantity is used to discard symmetry-breaking \(B_{mn}\) harmonics.
ntheta – Number of grid points in the poloidal angle for evaluating geometric quantities in the Redl formulae.
plot – Make a plot of many of the quantities computed.
- __call__()
Evaluate the geometric quantities needed for the Redl bootstrap current formula.
- class simsopt.mhd.RedlGeomVmec(vmec, surfaces=None, ntheta=64, nphi=65, plot=False)
Bases:
Optimizable
This class evaluates geometry data needed to evaluate the Redl bootstrap current formula from a vmec configuration, such as the effective fraction of trapped particles. The advantage of this class over
RedlGeomBoozer
is that no transformation to Boozer coordinates is involved in this method. However, the approach here may over-estimateepsilon
.- Parameters
vmec – An instance of
simsopt.mhd.vmec.Vmec
.surfaces – A 1D array of values of s (normalized toroidal flux) on which to compute the geometric quantities. If
None
, the half grid points from the VMEC solution will be used.ntheta – Number of grid points in the poloidal angle for evaluating geometric quantities in the Redl formulae.
nphi – Number of grid points in the toroidal angle for evaluating geometric quantities in the Redl formulae.
plot – Whether to make a plot of many of the quantities computed.
- __call__()
Evaluate the geometric quantities needed for the Redl bootstrap current formula.
- class simsopt.mhd.Residue(spec, pp, qq, vol=1, theta=0, s_guess=None, s_min=- 1.0, s_max=1.0, rtol=1e-09)
Bases:
Optimizable
Greene’s residue, evaluated from a Spec equilibrum
- Parameters
spec – a Spec object
pp – Numerator and denominator for the resonant iota = pp / qq
qq – Numerator and denominator for the resonant iota = pp / qq
vol – Index of the Spec volume to consider
theta – Spec’s theta coordinate at the periodic field line
s_guess – Guess for the value of Spec’s s coordinate at the periodic field line
s_min – bounds on s for the search
s_max – bounds on s for the search
rtol – the relative tolerance of the integrator
- J()
Run Spec if needed, find the periodic field line, and return the residue
- recompute_bell(parent=None)
Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.
This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.
Need to be implemented by classes that provide a dof_setter for external handling of DOFs.
- class simsopt.mhd.Spec(filename: Optional[str] = None, mpi: Optional[MpiPartition] = None, verbose: bool = True, keep_all_files: bool = False)
Bases:
Optimizable
This class represents the SPEC equilibrium code.
Philosophy regarding mpol and ntor: The Spec object keeps track of mpol and ntor values that are independent of those for the boundary Surface object. If the Surface object has different mpol/ntor values, the Surface’s rbc/zbs arrays are first copied, then truncated or expanded to fit the mpol/ntor values of the Spec object to before Spec is run. Therefore, you may sometimes need to manually change the mpol and ntor values for the Spec object.
The default behavior is that all output files will be deleted except for the first and most recent iteration on worker group 0. If you wish to keep all the output files, you can set
keep_all_files = True
. If you want to save the output files for a certain intermediate iteration, you can set thefiles_to_delete
attribute to[]
after that run of SPEC.- Parameters
filename – SPEC input file to use for initialization. It should end in
.sp
. Or, if None, default values will be used.mpi – A
simsopt.util.mpi.MpiPartition
instance, from which the worker groups will be used for SPEC calculations. IfNone
, each MPI process will run SPEC independently.verbose – Whether to print SPEC output to stdout.
keep_all_files – If
False
, all output files will be deleted except for the first and most recent ones from worker group 0. IfTrue
, all output files will be kept.
- property boundary
- get_dofs()
- init(filename: str)
Initialize SPEC fortran state from an input file.
- Parameters
filename – Name of the file to load. It should end in
.sp
.
- iota()
Return the rotational transform in the middle of the volume.
- 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.
- run()
Run SPEC, if needed.
- set_dofs(x)
- volume()
Return the volume inside the boundary flux surface.
- class simsopt.mhd.VirtualCasing
Bases:
object
Use the virtual casing principle to compute the contribution to the total magnetic field due to current outside a bounded surface.
Usually, an instance of this class is created using the
from_vmec()
class method, which also drives the computationally demanding part of the calculation (solving the integral equation). In the future, interfaces to other equilibrium codes may be added.In the standard 2-stage approach to stellarator optimization, the virtual casing calculation is run once, at the end of stage 1 (optimizing the plasma shape), with the result provided as input to stage 2 (optimizing the coil shapes). In this case, you can use the
save()
function or thefilename
argument offrom_vmec()
to save the results of the virtual casing calculation. These saved results can then be loaded in later using theload()
class method, when needed for solving the stage-2 problem.In order to compute the external field accurately, one requires fairly high grid resolution. For the stage-2 problem however, a lower resolution is often sufficient. To deal with this, we consider two grids, one denoted by the prefix
src_
and the other one denoted bytrgt_
.src_nphi
andsrc_ntheta
refer to the resolution of the grid for the input data, i.e. the total magnetic field and shape of the surface.trgt_nphi
andtrgt_ntheta
refer to the resolution of the grid where the external field is computed, i.e. the output of the virtual casing calculation that is provided as input to the stage 2 problem).To set the grid resolutions
src_nphi
andsrc_ntheta
, it can be convenient to use the functionsimsopt.geo.surface.best_nphi_over_ntheta()
.An instance of this class has the following attributes. For all vector quantites, Cartesian coordinates are used, corresponding to array dimensions of size 3:
src_nphi
: The number of grid points in the toroidal angle \(\phi\), for a half field period (or a full field period if use_stellsym=False)src_ntheta
: The number of grid points in the poloidal angle \(\theta\).src_phi
: An array of size(src_nphi,)
with the grid points of \(\phi\).src_theta
: An array of size(src_ntheta,)
with the grid points of \(\theta\).trgt_nphi
: The number of grid points in the toroidal angle \(\phi\), for a half field period (or a full field period if use_stellsym=False)trgt_ntheta
: The number of grid points in the poloidal angle \(\theta\).trgt_phi
: An array of size(trgt_nphi,)
with the grid points of \(\phi\).trgt_theta
: An array of size(trgt_ntheta,)
with the grid points of \(\theta\).gamma
: An array of size(src_nphi, src_ntheta, 3)
with the position vector on the surface.B_total
: An array of size(src_nphi, src_ntheta, 3)
with the total magnetic field vector on the surface.unit_normal
: An array of size(trgt_nphi, trgt_ntheta, 3)
with the unit normal vector on the surface.B_external
: An array of size(trgt_nphi, trgt_ntheta, 3)
with the contribution to the magnetic field due to current outside the surface.B_external_normal
: An array of size(trgt_nphi, trgt_ntheta)
with the contribution to the magnetic field due to current outside the surface, taking just the component normal to the surface.
The \(\phi\) and \(\theta\) grids for these data are both uniformly spaced, and are the same as for
Surface
classes withrange="half period"
orrange="full period"
, for the case of stellarator-symmetry or non-stellarator-symmetry respectively. (See the description of therange
parameter in the documentation on Surfaces.) For the usual case of stellarator symmetry, all the virtual casing data are given on half a field period. There is no grid point at \(\phi=0\), rather the grid is shifted in \(\phi\) by half the grid spacing. Thus, thesrc_phi
grid isnp.linspace(1 / (2 * nfp * src_nphi), (src_nphi - 0.5) / (src_nphi * nfp), src_nphi)
(recalling the simsopt convention that \(\phi\) and \(\theta\) have period 1, not \(2\pi\)). For a non-stellarator-symmetric calculation, thesrc_phi
grid isnp.linspace(0, 1 / nfp, src_nphi, endpoint=False)
. Thetrgt_phi
grid follows the same logic as thesrc_phi
grid. Note that for stellarator symmetry, ifsrc_nphi != trgt_nphi
, then the shift (i.e. first grid point) insrc_phi
andtrgt_phi
will be different. For both stellarator symmetry and non-stellarator-symmetry, thesrc_theta
grid isnp.linspace(0, 1, src_ntheta, endpoint=False)
, and thetrgt_theta
grid is the same but withtrgt_ntheta
.In particular,
B_external_normal
is given on the grid that would be naturally used for stage-2 coil optimization, so no resampling is required.- classmethod from_vmec(vmec, src_nphi, src_ntheta=None, trgt_nphi=None, trgt_ntheta=None, use_stellsym=True, digits=6, filename='auto')
Given a
Vmec
object, compute the contribution to the total magnetic field due to currents outside the plasma.This function requires the python
virtual_casing
package to be installed.The argument
src_nphi
refers to the number of points around a half field period if stellarator symmetry is exploited, or a full field period if not.To set the grid resolutions
src_nphi
andsrc_ntheta
, it can be convenient to use the functionsimsopt.geo.surface.best_nphi_over_ntheta()
. This is done automatically if you omit thesrc_ntheta
argument.For now, this routine only works for stellarator symmetry.
- Parameters
vmec – Either an instance of
simsopt.mhd.vmec.Vmec
, or the name of a Vmecinput.*
orwout*
file.src_nphi – Number of grid points toroidally for the input of the calculation.
src_ntheta – Number of grid points poloidally for the input of the calculation. If
None
, the number of grid points will be calculated automatically usingsimsopt.geo.surface.best_nphi_over_ntheta()
to minimize the grid anisotropy, given the specifiednphi
.trgt_nphi – Number of grid points toroidally for the output of the calculation. If unspecified,
src_nphi
will be used.trgt_ntheta – Number of grid points poloidally for the output of the calculation. If unspecified,
src_ntheta
will be used.use_stellsym – whether to exploit stellarator symmetry in the calculation.
digits – Approximate number of digits of precision for the calculation.
filename – If not
None
, the results of the virtual casing calculation will be saved in this file. For the default value of"auto"
, the filename will automatically be set to"vcasing_<extension>.nc"
where<extension>
is the string associated with Vmec input and output files, analogous to the Vmec output file"wout_<extension>.nc"
.
- classmethod load(filename)
Load in the results of a previous virtual casing calculation, previously saved in NetCDF format.
- Parameters
filename – Name of the file to load.
- plot(ax=None, show=True)
Plot
B_external_normal
, the component normal to the surface of the magnetic field generated by currents outside the surface. This routine requiresmatplotlib
.- Parameters
ax – The axis object on which to plot. This argument is useful when plotting multiple objects on the same axes. If equal to the default
None
, a new axis will be created.show – Whether to call matplotlib’s
show()
function.
- Returns
An axis which could be passed to a further call to matplotlib if desired.
- save(filename='vcasing.nc')
Save the results of a virtual casing calculation in a NetCDF file.
- Parameters
filename – Name of the file to create.
- class simsopt.mhd.Vmec(filename: Optional[str] = None, mpi: Optional[MpiPartition] = None, keep_all_files: bool = False, verbose: bool = True, ntheta=50, nphi=50)
Bases:
Optimizable
This class represents the VMEC equilibrium code.
You can initialize this class either from a VMEC
input.<extension>
file or from awout_<extension>.nc
output file. If neither is provided, a default input file is used. When this class is initialized from an input file, it is possible to modify the input parameters and run the VMEC code. When this class is initialized from awout
file, all the data from thewout
file is available in memory but the VMEC code cannot be re-run, since some of the input data (e.g. radial multigrid parameters) is not available in the wout file.The input parameters to VMEC are all accessible as attributes of the
indata
attribute. For example, ifvmec
is an instance ofVmec
, then you can read or write the input resolution parameters usingvmec.indata.mpol
,vmec.indata.ntor
,vmec.indata.ns_array
, etc. However, the boundary surface is different:rbc
,rbs
,zbc
, andzbs
from theindata
attribute are always ignored, and these arrays are instead taken from the simsopt surface object associated to theboundary
attribute. Ifboundary
is a surface based on some other representation than VMEC’s Fourier representation, the surface will automatically be converted to VMEC’s representation (SurfaceRZFourier
) before each run of VMEC. You can replaceboundary
with a new surface object, of any type that implements the conversion functionto_RZFourier()
.VMEC is run either when the
run()
function is called, or when any of the output functions likeaspect()
oriota_axis()
are called.A caching mechanism is implemented, using the attribute
need_to_run_code
. Whenever VMEC is run, or if the class is initialized from awout
file, this attribute is set toFalse
. Subsequent calls torun()
or output functions likeaspect()
will not actually run VMEC again, untilneed_to_run_code
is changed toTrue
. The attributeneed_to_run_code
is automatically set toTrue
whenever the state vector.x
is changed, and when dofs of theboundary
are changed. However,need_to_run_code
is not automatically set toTrue
when entries ofindata
are modified.Once VMEC has run at least once, or if the class is initialized from a
wout
file, all of the quantities in thewout
output file are available as attributes of thewout
attribute. For example, ifvmec
is an instance ofVmec
, then the flux surface shapes can be obtained fromvmec.wout.rmnc
andvmec.wout.zmns
.Since the underlying fortran implementation of VMEC uses global module variables, it is not possible to have more than one python Vmec object with different parameters; changing the parameters of one would change the parameters of the other.
An instance of this class owns just a few optimizable degrees of freedom, particularly
phiedge
andcurtor
. The optimizable degrees of freedom associated with the boundary surface are owned by that surface object.To run VMEC, two input profiles must be specified: pressure and either iota or toroidal current. Each of these profiles can be specified in several ways. One way is to specify the profile in the input file used to initialize the
Vmec
object. For instance, the pressure profile is determined by the variablespmass_type
,am
,am_aux_s
, andam_aux_f
. You can also modify these variables from python via theindata
attribute, e.g.vmec.indata.am = [1.0e5, -1.0e5]
. Another option is to assign asimsopt.mhd.profiles.Profile
object to the attributespressure_profile
,current_profile
, oriota_profile
. This approach allows for the profiles to be optimized, and it allows you to use profile shapes defined in python that are not available in the fortran VMEC code. To explain this approach we focus here on the pressure profile; the iota and current profiles are analogous. If thepressure_profile
attribute of aVmec
object isNone
(the default), then a simsoptProfile
object is not used, and instead the settings fromVmec.indata
(initialized from the input file) are used. If aProfile
object is assigned to thepressure_profile
attribute, then an edge in the dependency graph is introduced, so theVmec
object then depends on the dofs of theProfile
object. Whenever VMEC is run, the simsoptProfile
is converted to either a polynomial (power series) or cubic spline in the normalized toroidal flux \(s\), depending on whetherindata.pmass_type
is"power_series"
or"cubic_spline"
. (The current profile is different in that either"cubic_spline_ip"
or"cubic_spline_i"
is specified instead of"cubic_spline"
.) The number of terms in the power series or number of spline nodes is determined by the attributesn_pressure
,n_current
, andn_iota
. If a cubic spline is used, the spline nodes are uniformly spaced from \(s=0\) to 1. Note that the choice of whether a polynomial or spline is used for the VMEC calculation is independent of the subclass ofProfile
used. Also, whether the iota or current profile is used is always determined by theindata.ncurr
attribute: 0 for iota, 1 for current. Example:from sismopt.mhd.profiles import ProfilePolynomial, ProfileSpline, ProfilePressure, ProfileScaled from simsopt.util.constants import ELEMENTARY_CHARGE ne = ProfilePolynomial(1.0e20 * np.array([1, 0, 0, 0, -0.9])) Te = ProfilePolynomial(8.0e3 * np.array([1, -0.9])) Ti = ProfileSpline([0, 0.5, 0.8, 1], 7.0e3 * np.array([1, 0.9, 0.8, 0.1])) ni = ne pressure = ProfilePressure(ne, Te, ni, Ti) # p = ne * Te + ni * Ti pressure_Pa = ProfileScaled(pressure, ELEMENTARY_CHARGE) # Te and Ti profiles were in eV, so convert to SI here. vmec = Vmec(filename) vmec.pressure_profile = pressure_Pa vmec.indata.pmass_type = "cubic_spline" vmec.n_pressure = 8 # Use 8 spline nodes
When VMEC is run multiple times, the default behavior is that all
wout
output files will be deleted except for the first and most recent iteration on worker group 0. If you wish to keep all thewout
files, you can setkeep_all_files = True
. If you want to save thewout
file for a certain intermediate iteration, you can set thefiles_to_delete
attribute to[]
after that run of VMEC.- Parameters
filename – Name of a VMEC
input.<extension>
file orwout_<extension>.nc
output file to use for loading the initial parameters. IfNone
, default parameters will be used.mpi – A
simsopt.util.mpi.MpiPartition
instance, from which the worker groups will be used for VMEC calculations. IfNone
, each MPI process will run VMEC independently.keep_all_files – If
False
, allwout
output files will be deleted except for the first and most recent ones from worker group 0. IfTrue
, allwout
files will be kept.verbose – Whether to print to stdout when running vmec.
- iter
Number of times VMEC has run.
- s_full_grid
The “full” grid in the radial coordinate s (normalized toroidal flux), including points at s=0 and s=1. Used for the output arrays and
zmns
.
- s_half_grid
The “half” grid in the radial coordinate s, used for
bmnc
,lmns
, and other output arrays. In contrast to wout files, this array has only ns-1 entries, so there is no leading 0.
- ds
The spacing between grid points for the radial coordinate s.
- __repr__()
Print the object in an informative way.
- aspect()
Return the plasma aspect ratio.
- property boundary
- property current_profile
- external_current()
Return the total electric current associated with external currents, i.e. the current through the “doughnut hole”. This number is useful for coil optimization, to know what the sum of the coil currents must be.
- Returns
float with the total external electric current in Amperes.
- get_dofs()
- get_input()
Generate a VMEC input file. The result will be returned as a string. To save a file, see the
write_input()
function.
- get_max_mn()
Look through the rbc and zbs data in fortran to determine the largest m and n for which rbc or zbs is nonzero.
- iota_axis()
Return the rotational transform on axis
- iota_edge()
Return the rotational transform at the boundary
- property iota_profile
- load_wout()
Read in the most recent
wout
file created, and store all the data in awout
attribute of this Vmec object.
- mean_iota()
Return the mean rotational transform. The average is taken over the normalized toroidal flux s.
- mean_shear()
Return an average magnetic shear, d(iota)/ds, where s is the normalized toroidal flux. This is computed by fitting the rotational transform to a linear (plus constant) function in s. The slope of this fit function is returned.
- property pressure_profile
- recompute_bell(parent=None)
Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.
This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.
Need to be implemented by classes that provide a dof_setter for external handling of DOFs.
- return_fn_map: Dict[str, Callable] = {'aspect': <function Vmec.aspect>, 'iota_axis': <function Vmec.iota_axis>, 'iota_edge': <function Vmec.iota_edge>, 'mean_iota': <function Vmec.mean_iota>, 'mean_shear': <function Vmec.mean_shear>, 'vacuum_well': <function Vmec.vacuum_well>, 'volume': <function Vmec.volume>}
- run()
Run VMEC, if
need_to_run_code
isTrue
.
- set_dofs(x)
- set_indata()
Transfer data from simsopt objects to Vmec’s fortran module data. Presently, this function sets the boundary shape and magnetic axis shape. In the future, the input profiles will be set here as well. This data transfer is performed before writing a Vmec input file or running Vmec. The boundary surface object converted to
SurfaceRZFourier
is returned.
- set_profile(longname, shortname, letter)
This function is used to set the pressure, current, and/or iota profiles.
- update_mpi(new_mpi)
Replace the
MpiPartition
with a new one.- Parameters
new_mpi – A new
simsopt.util.mpi.MpiPartition
object.
- vacuum_well()
Compute a single number W that summarizes the vacuum magnetic well, given by the formula
W = (dV/ds(s=0) - dV/ds(s=1)) / (dV/ds(s=0)
where dVds is the derivative of the flux surface volume with respect to the radial coordinate s. Positive values of W are favorable for stability to interchange modes. This formula for W is motivated by the fact that
d^2 V / d s^2 < 0
is favorable for stability. Integrating over s from 0 to 1 and normalizing gives the above formula for W. Notice that W is dimensionless, and it scales as the square of the minor radius. To compute dV/ds, we use
dV/ds = 4 * pi**2 * abs(sqrt(g)_{0,0})
where sqrt(g) is the Jacobian of (s, theta, phi) coordinates, computed by VMEC in the gmnc array, and _{0,0} indicates the m=n=0 Fourier component. Since gmnc is reported by VMEC on the half mesh, we extrapolate by half of a radial grid point to s = 0 and 1.
- volume()
Return the volume inside the VMEC last closed flux surface.
- write_input(filename)
Write a VMEC input file. To just get the result as a string without saving a file, see the
get_input()
function.- Parameters
filename – Name of the file to write. Selected MPI processes can pass
None
if you wish for these processes to not write a file.
- class simsopt.mhd.VmecRedlBootstrapMismatch(geom, ne, Te, Ti, Zeff, helicity_n, logfile=None)
Bases:
Optimizable
This class is used to obtain quasi-axisymmetric or quasi-helically symmetric VMEC configurations with self-consistent bootstrap current. This class represents the objective function
\[f = \frac{\int ds \left[\left<\vec{J}\cdot\vec{B}\right>_{vmec} - \left<\vec{J}\cdot\vec{B}\right>_{Redl} \right]^2} {\int ds \left[\left<\vec{J}\cdot\vec{B}\right>_{vmec} + \left<\vec{J}\cdot\vec{B}\right>_{Redl} \right]^2}\]where \(\left<\vec{J}\cdot\vec{B}\right>_{vmec}\) is the bootstrap current profile in a VMEC equilibrium, and \(\left<\vec{J}\cdot\vec{B}\right>_{Redl}\) is the bootstrap current profile computed from the fit formulae in Redl et al, Physics of Plasmas 28, 022502 (2021).
- Parameters
geom – An instance of either
RedlGeomVmec
orRedlGeomBoozer
.ne – A
Profile
object representing the electron density profile.Te – A
Profile
object representing the electron temperature profile.Ti – A
Profile
object representing the ion temperature profile.Zeff – A
Profile
object representing the \(Z_{eff}\) profile. A single number can also be provided, in which case a constant \(Z_{eff}\) profile will be used.helicity_n – 0 for quasi-axisymmetry, or +/- 1 for quasi-helical symmetry.
- J()
Return the scalar objective function, given by the sum of the squares of the residuals.
- residuals()
This function returns a 1d array of residuals, useful for representing the objective function as a nonlinear least-squares problem. This is the function handle to use with a
LeastSquaresProblem
.Specifically, this function returns
\[R_j = \frac{\left<\vec{J}\cdot\vec{B}\right>_{vmec}(s_j) - \left<\vec{J}\cdot\vec{B}\right>_{Redl}(s_j)} {\sqrt{\sum_{k=1}^N \left[\left<\vec{J}\cdot\vec{B}\right>_{vmec}(s_k) + \left<\vec{J}\cdot\vec{B}\right>_{Redl}(s_k) \right]^2}}\]where \(j\) and \(k\) range over the surfaces for the supplied
geom
object (typically the half-grid points for the VMEC configuration), \(j, k \in \{1, 2, \ldots, N\}\) and \(N\) is the number of surfaces for the suppliedgeom
object. This corresponds to approximating the \(\int ds\) integrals in the objective function with Riemann integration. The vector of residuals returned has length \(N\).The sum of the squares of these residuals equals the objective function. The total scalar objective is approximately independent of the number of surfaces.
- class simsopt.mhd.WellWeighted(vmec, weight_function1, weight_function2, adjoint_epsilon=0.1)
Bases:
Optimizable
WellWeighted is a class that computes a measure of magnetic well for a vmec equilibrium. The magnetic well measure is
\[J = \frac{ \int ds \, V'(s) [w_1(s) - w_2(s)]} { \int ds \, V'(s) [w_1(s) + w_2(s)]},\]where \(w_1(s)\) and \(w_2(s)\) correspond to the arguments
weight_function1
andweight_function2
, and \(V(s)\) is the volume enclosed by the flux surface with normalized toroidal flux \(s\). Typically, \(w_1\) would be peaked on the edge while \(w_2\) would be peaked on the axis, such that \(J < 0\) corresonds to \(V''(s) < 0\), which is favorable for stability.This class also provides calculations of the derivatives of \(J\) using an adjoint method.
- Parameters
vmec – instance of Vmec
weight_function1 – function handle which takes a single argument, s, the normalized toroidal flux
weight_function2 – function handle which takes a single argument, s, the normalized toroidal flux
adjoint_epsilon – sets the amplitude of the toroidal current perturbation required for the adjoint solve.
- J()
Computes the quantity \(J\) described in the class definition.
- dJ()
Computes derivatives of \(J\) with respect to surface parameters using an adjoint method.
- shape_gradient()
Computes the shape gradient of the quantity \(J\) described in the class definition. For a perturbation to the surface \(\delta \vec{x}\), the resulting perturbation to the objective function is
\[\delta J(\delta \vec{x}) = \int d^2 x \, G \delta \vec{x} \cdot \vec{n}\]where the integral is over the VMEC boundary surface, \(G\) is the shape gradient, and \(\vec{n}\) is the unit normal.
- Returns
2d array of size (numquadpoints_phi,numquadpoints_theta)
- Return type
\(G\)
- simsopt.mhd.compute_trapped_fraction(modB, sqrtg)
Compute the effective fraction of trapped particles, which enters several formulae for neoclassical transport, as well as several quantities that go into its calculation. The input data can be provided on a uniform grid of arbitrary toroidal and poloidal angles that need not be straight-field-line angles.
The trapped fraction
f_t
has a standard definition in neoclassical theory:\[f_t = 1 - \frac{3}{4} \left< B^2 \right> \int_0^{1/Bmax} \frac{\lambda\; d\lambda}{\left< \sqrt{1 - \lambda B} \right>}\]where \(\left< \ldots \right>\) is a flux surface average.
The effective inverse aspect ratio epsilon is defined by
\[\frac{Bmax}{Bmin} = \frac{1 + \epsilon}{1 - \epsilon}\]This definition is motivated by the fact that this formula would be true in the case of circular cross-section surfaces in axisymmetry with \(B \propto 1/R\) and \(R = (1 + \epsilon \cos\theta) R_0\).
- Parameters
modB – 2D array of size (ntheta, ns) or 3D array of size (ntheta, nphi, ns) with \(|B|\) on the grid points.
sqrtg – 2D array of size (ntheta, ns) or 3D array of size (ntheta, nphi, ns) with the Jacobian \(1/(\nabla s \times\nabla\theta\cdot\nabla\phi)\) on the grid points.
- Returns
Tuple containing
Bmin: A 1D array, with the minimum of \(|B|\) on each surface.
Bmax: A 1D array, with the maximum of \(|B|\) on each surface.
epsilon: A 1D array, with the effective inverse aspect ratio on each surface.
fsa_B2: A 1D array with \(\left<B^2\right>\) on each surface, where \(\left< \ldots \right>\) denotes a flux surface average.
fsa_1overB: A 1D array with \(\left<1/B\right>\) on each surface, where \(\left< \ldots \right>\) denotes a flux surface average.
f_t: A 1D array, with the effective trapped fraction on each surface.
- simsopt.mhd.j_dot_B_Redl(ne, Te, Ti, Zeff, helicity_n=None, s=None, G=None, R=None, iota=None, epsilon=None, f_t=None, psi_edge=None, nfp=None, geom=None, plot=False)
Compute the bootstrap current (specifically \(\left<\vec{J}\cdot\vec{B}\right>\)) using the formulae in Redl et al, Physics of Plasmas 28, 022502 (2021).
The profiles of ne, Te, Ti, and Zeff should all be instances of subclasses of
simsopt.mhd.profiles.Profile
, i.e. they should have__call__()
anddfds()
functions. IfZeff == None
, a constant 1 is assumed. IfZeff
is a float, a constant profile will be assumed.ne
should have units of 1/m^3.Ti
andTe
should have units of eV.Geometric data can be specified in one of two ways. In the first approach, the arguments
s
,G
,R
,iota
,epsilon
,f_t
,psi_edge
, andnfp
are specified, while the argumentgeom
is not. In the second approach, the argumentgeom
is set to an instance of eitherRedlGeomVmec
orRedlGeomBoozer
, and this object will be used to set all the other geometric quantities. In this case, the argumentss
,G
,R
,iota
,epsilon
,f_t
,psi_edge
, andnfp
should not be specified.The input variable
s
is a 1D array of values of normalized toroidal flux. The input arraysG
,R
,iota
,epsilon
, andf_t
, should be 1d arrays evaluated on this sames
grid. The bootstrap current \(\left<\vec{J}\cdot\vec{B}\right>\) will be computed on this same set of flux surfaces.If you provide a
RedlGeomBoozer
object forgeom
, then it is not necessary to specify the argumenthelicity_n
here, in which casehelicity_n
will be taken fromgeom
.- Parameters
ne – A
Profile
object with the electron density profile.Te – A
Profile
object with the electron temperature profile.Ti – A
Profile
object with the ion temperature profile.Zeff – A
Profile
object with the profile of the average impurity charge \(Z_{eff}\). Or, a single number can be provided if this profile is constant. Or, ifNone
, Zeff = 1 will be used.helicity_n – 0 for quasi-axisymmetry, or +/- 1 for quasi-helical symmetry. This quantity is used to apply the quasisymmetry isomorphism to map the collisionality and bootstrap current from the tokamak expressions to quasi-helical symmetry.
s – A 1D array of values of normalized toroidal flux.
G – A 1D array with the flux function multiplying \(\nabla\varphi\) in the Boozer covariant representation, equivalent to \(R B_{toroidal}\) in axisymmetry.
R – A 1D array with the effective major radius to use when evaluating the collisionality in the Sauter/Redl formulae.
iota – A 1D array with the rotational transform.
epsilon – A 1D array with the effective inverse aspect ratio to use for evaluating the collisionality in the Sauter/Redl formulae.
f_t – A 1D array with the effective trapped fraction.
psi_edge – The toroidal flux (in Webers) divided by (2pi) at the boundary s=1
nfp – The number of field periods. Irrelevant for axisymmetry or quasi-axisymmetry; matters only if
helicity_n
is not 0.geom – Optional. An instance of either
RedlGeomVmec
orRedlGeomBoozer
.plot – Whether to make a plot of many of the quantities computed.
- Returns
Tuple containing
jdotB: A 1D array containing the bootstrap current \(\left<\vec{J}\cdot\vec{B}\right>\) on the specified flux surfaces.
details: An object holding intermediate quantities from the computation (e.g. L31, L32, alpha) as attributes