simsopt.geo package
Submodules
simsopt.geo.boozersurface module
- class simsopt.geo.boozersurface.BoozerSurface(biotsavart, surface, label, targetlabel)
Bases:
object
BoozerSurface and its associated methods can be used to compute the Boozer angles on a surface. It takes a Surface representation (e.g. SurfaceXYZFourier, or SurfaceXYZTensorFourier), a magnetic field evaluator, surface label evaluator, and a target surface label.
The Boozer angles are computed by solving a constrained least squares problem. The least squares objective is given by \(J(x) = \frac{1}{2} \mathbf r^T(x) \mathbf r(x)\), where \(\mathbf r\) is a vector of residuals computed by
boozer_surface_residual
(seesurfaceobjectives.py
), and some constraints. This objective is zero when the surface corresponds to a magnetic surface of the field and \((\phi,\theta)\) that parametrize the surface correspond to Boozer angles, and the constraints are satisfied.The surface label can be area, volume, or toroidal flux. The label on the computed surface will be equal or close to the user-provided
targetlabel
, depending on how the label constraint is imposed. This constrained least squares problem can be solved by scalarizing and adding the constraint as an additional penalty term to the objective. This is done inminimize_boozer_penalty_constraints_LBFGS
minimize_boozer_penalty_constraints_newton
minimize_boozer_penalty_constraints_ls
where LBFGS, Newton, or
scipy.optimize.least_squares
optimizers are used, respectively. Alternatively, the exactly constrained least squares optimization problem can be solved. This is done inminimize_boozer_exact_constraints_newton
where Newton is used to solve the first order necessary conditions for optimality.
- boozer_exact_constraints(xl, derivatives=0, optimize_G=True)
This function returns the optimality conditions corresponding to the minimization problem
\[\text{min}_x ~J(x)\]subject to
\[\begin{split}l - l_0 &= 0 \\ z(\varphi=0,\theta=0) - 0 &= 0\end{split}\]where \(l\) is the surface label and \(l_0\) is the target surface label, \(J(x) = \frac{1}{2}\mathbf r(x)^T \mathbf r(x)\), and \(\mathbf r(x)\) contains the Boozer residuals at quadrature points \(1,\dots,n\). We can also optionally return the first derivatives of these optimality conditions.
- boozer_penalty_constraints(x, derivatives=0, constraint_weight=1.0, scalarize=True, optimize_G=False)
Define the residual
\[\mathbf r(x) = [r_1(x),...,r_n(x), \sqrt{w_c} (l-l_0), \sqrt{w_c} (z(\varphi=0, \theta=0) - 0)]\]where \(w_c\) is the constraint weight, \(r_i\) are the Boozer residuals at quadrature points \(1,\dots,n\), \(l\) is the surface label, and \(l_0\) is the target surface label.
For
scalarized=False
, this function returns \(\mathbf r(x)\) and optionally the Jacobian of \(\mathbf r(x)\).for
scalarized=True
, this function returns\[J(x) = \frac{1}{2}\mathbf r(x)^T \mathbf r(x),\]i.e. the least squares residual and optionally the gradient and the Hessian of \(J(x)\).
- minimize_boozer_exact_constraints_newton(tol=1e-12, maxiter=10, iota=0.0, G=None, lm=[0.0, 0.0])
This function solves the constrained optimization problem
\[\text{min}_x ~ J(x)\]subject to
\[\begin{split}l - l_0 &= 0 \\ z(\varphi=0,\theta=0) - 0 &= 0\end{split}\]using Lagrange multipliers and Newton’s method. In the above, \(J(x) = \frac{1}{2}\mathbf r(x)^T \mathbf r(x)\), and \(\mathbf r(x)\) contains the Boozer residuals at quadrature points \(1,\dots,n\).
The final constraint is not necessary for stellarator symmetric surfaces as it is automatically satisfied by the stellarator symmetric surface parametrization.
- minimize_boozer_penalty_constraints_LBFGS(tol=0.001, maxiter=1000, constraint_weight=1.0, iota=0.0, G=None)
This function tries to find the surface that approximately solves
\[\text{min}_x ~J(x) + \frac{1}{2} w_c (l - l_0)^2 + \frac{1}{2} w_c (z(\varphi=0, \theta=0) - 0)^2\]where \(J(x) = \frac{1}{2}\mathbf r(x)^T \mathbf r(x)\), and \(\mathbf r(x)\) contains the Boozer residuals at quadrature points \(1,\dots,n\). This is done using LBFGS.
- minimize_boozer_penalty_constraints_ls(tol=1e-12, maxiter=10, constraint_weight=1.0, iota=0.0, G=None, method='lm')
This function does the same as
minimize_boozer_penalty_constraints_LBFGS
, but instead of LBFGS it uses a nonlinear least squares algorithm whenmethod='lm'
. Options for the method are the same as forscipy.optimize.least_squares
. Ifmethod='manual'
, then a damped Gauss-Newton method is used.
- minimize_boozer_penalty_constraints_newton(tol=1e-12, maxiter=10, constraint_weight=1.0, iota=0.0, G=None, stab=0.0)
This function does the same as
minimize_boozer_penalty_constraints_LBFGS
, but instead of LBFGS it uses Newton’s method.
- solve_residual_equation_exactly_newton(tol=1e-10, maxiter=10, iota=0.0, G=None)
This function solves the Boozer Surface residual equation exactly. For this to work, we need the right balance of quadrature points, degrees of freedom and constraints. For this reason, this only works for surfaces of type SurfaceXYZTensorFourier right now.
Given ntor, mpol, nfp and stellsym, the surface is expected to be created in the following way:
phis = np.linspace(0, 1/nfp, 2*ntor+1, endpoint=False) thetas = np.linspace(0, 1, 2*mpol+1, endpoint=False) s = SurfaceXYZTensorFourier(
mpol=mpol, ntor=ntor, stellsym=stellsym, nfp=nfp, quadpoints_phi=phis, quadpoints_theta=thetas)
- Or the following two are also possible in the stellsym case
phis = np.linspace(0, 1/nfp, 2*ntor+1, endpoint=False) thetas = np.linspace(0, 0.5, mpol+1, endpoint=False)
- or
phis = np.linspace(0, 1/(2*nfp), ntor+1, endpoint=False) thetas = np.linspace(0, 1, 2*mpol+1, endpoint=False)
and then
- s = SurfaceXYZTensorFourier(
mpol=mpol, ntor=ntor, stellsym=stellsym, nfp=nfp, quadpoints_phi=phis, quadpoints_theta=thetas)
For the stellsym case, there is some redundancy between dofs. This is taken care of inside this function.
- Not stellsym:
The surface has (2*ntor+1)*(2*mpol+1) many quadrature points and 3*(2*ntor+1)*(2*mpol+1) many dofs. Equations:
Boozer residual in x, y, and z at all quadrature points
z(0, 0) = 0
label constraint (e.g. volume or flux)
- Unknowns:
Surface dofs
iota
G
So we end up having 3*(2*ntor+1)*(2*mpol+1) + 2 equations and the same number of unknowns.
- Stellsym:
In this case we have
- D = (ntor+1)*(mpol+1)+ ntor*mpol + 2*(ntor+1)*mpol + 2*ntor*(mpol+1)
= 6*ntor*mpol + 3*ntor + 3*mpol + 1
many dofs in the surface. After calling surface.get_stellsym_mask() we have kicked out
2*ntor*mpol + ntor + mpol
quadrature points, i.e. we have
2*ntor*mpol + ntor + mpol + 1
quadrature points remaining. In addition we know that the x coordinate of the residual at phi=0=theta is also always satisfied. In total this leaves us with
3*(2*ntor*mpol + ntor + mpol) + 2 equations for the boozer residual. 1 equation for the label
which is the same as the number of surface dofs + 2 extra unknowns given by iota and G.
simsopt.geo.config module
simsopt.geo.curve module
- class simsopt.geo.curve.Curve(**kwargs)
Bases:
simsopt._core.graph_optimizable.Optimizable
Curve is a base class for various representations of curves in SIMSOPT using the graph based Optimizable framework with external handling of DOFS as well.
- dfrenet_frame_by_dcoeff()
This function returns the derivative of the curve’s Frenet frame,
\[\left(\frac{\partial \mathbf{t}}{\partial \mathbf{c}}, \frac{\partial \mathbf{n}}{\partial \mathbf{c}}, \frac{\partial \mathbf{b}}{\partial \mathbf{c}}\right),\]with respect to the curve dofs, where \((\mathbf t, \mathbf n, \mathbf b)\) correspond to the Frenet frame, and \(\mathbf c\) are the curve dofs.
- dgamma_by_dcoeff_vjp(v)
- dgammadash_by_dcoeff_vjp(v)
- dgammadashdash_by_dcoeff_vjp(v)
- dgammadashdashdash_by_dcoeff_vjp(v)
- dincremental_arclength_by_dcoeff_vjp(v)
This function returns the vector Jacobian product
\[v^T \frac{\partial \|\Gamma'\|}{\partial \mathbf{c}}\]where \(\|\Gamma'\|\) is the incremental arclength, \(\Gamma'\) is the tangent to the curve and \(\mathbf{c}\) are the curve dofs.
- dkappa_by_dcoeff_impl(dkappa_by_dcoeff)
This function computes the derivative of the curvature with respect to the curve coefficients.
\[\frac{\partial \kappa}{\partial \mathbf c}\]where \(\mathbf c\) are the curve dofs, \(\kappa\) is the curvature.
- dkappa_by_dcoeff_vjp(v)
This function returns the vector Jacobian product
\[v^T \frac{\partial \kappa}{\partial \mathbf{c}}\]where \(\mathbf c\) are the curve dofs and \(\kappa\) is the curvature.
- dkappadash_by_dcoeff()
This function returns
\[\frac{\partial \kappa'(\phi)}{\partial \mathbf{c}}.\]where \(\mathbf{c}\) are the curve dofs, and \(\kappa\) is the curvature.
- dtorsion_by_dcoeff_impl(dtorsion_by_dcoeff)
This function returns the derivative of torsion with respect to the curve dofs.
\[\frac{\partial \tau}{\partial \mathbf c}\]where \(\mathbf c\) are the curve dofs, and \(\tau\) is the torsion.
- dtorsion_by_dcoeff_vjp(v)
This function returns the vector Jacobian product
\[v^T \frac{\partial \tau}{\partial \mathbf{c}}\]where \(\mathbf c\) are the curve dofs, and \(\tau\) is the torsion.
- frenet_frame()
This function returns the Frenet frame, \((\mathbf{t}, \mathbf{n}, \mathbf{b})\), associated to the curve.
- kappa_impl(kappa)
This function implements the curvature, \(\kappa(\varphi)\).
- kappadash()
This function returns \(\kappa'(\phi)\), where \(\kappa\) is the curvature.
- plot(engine='matplotlib', ax=None, show=True, plot_derivative=False, close=False, axis_equal=True, **kwargs)
Plot the curve in 3D using
matplotlib.pyplot
,mayavi
, orplotly
.- Parameters
engine – The graphics engine to use. Available settings are
"matplotlib"
,"mayavi"
, and"plotly"
.ax – The axis object on which to plot. This argument is useful when plotting multiple objects on the same axes. If equal to the default
None
, a new axis will be created.show – Whether to call the
show()
function of the graphics engine. Should be set toFalse
if more objects will be plotted on the same axes.plot_derivative – Whether to plot the tangent of the curve too. Not implemented for plotly.
close – Whether to connect the first and last point on the curve. Can lead to surprising results when only quadrature points on a part of the curve are considered, e.g. when exploting rotational symmetry.
axis_equal – For matplotlib, whether all three dimensions should be scaled equally.
kwargs – Any additional arguments to pass to the plotting function, like
color='r'
.
- Returns
An axis which could be passed to a further call to the graphics engine so multiple objects are shown together.
- recompute_bell(parent=None)
For derivative classes of Curve, all of which also subclass from C++ Curve class, call invalidate_cache which is implemented in C++ side.
- torsion_impl(torsion)
This function returns the torsion, \(\tau\), of a curve.
- class simsopt.geo.curve.JaxCurve(quadpoints, gamma_pure, **kwargs)
Bases:
simsoptpp.Curve
,simsopt.geo.curve.Curve
- dgamma_by_dcoeff_impl(dgamma_by_dcoeff)
This function returns
\[\frac{\partial \Gamma}{\partial \mathbf c}\]where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.
- dgamma_by_dcoeff_vjp_impl(v)
This function returns the vector Jacobian product
\[v^T \frac{\partial \Gamma}{\partial \mathbf c}\]where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.
- dgammadash_by_dcoeff_impl(dgammadash_by_dcoeff)
This function returns
\[\frac{\partial \Gamma'}{\partial \mathbf c}\]where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.
- dgammadash_by_dcoeff_vjp_impl(v)
This function returns
\[\mathbf v^T \frac{\partial \Gamma'}{\partial \mathbf c}\]where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.
- dgammadashdash_by_dcoeff_impl(dgammadashdash_by_dcoeff)
This function returns
\[\frac{\partial \Gamma''}{\partial \mathbf c}\]where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.
- dgammadashdash_by_dcoeff_vjp_impl(v)
This function returns the vector Jacobian product
\[v^T \frac{\partial \Gamma''}{\partial \mathbf c}\]where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.
- dgammadashdashdash_by_dcoeff_impl(dgammadashdashdash_by_dcoeff)
This function returns
\[\frac{\partial \Gamma'''}{\partial \mathbf c}\]where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.
- dgammadashdashdash_by_dcoeff_vjp_impl(v)
This function returns the vector Jacobian product
\[v^T \frac{\partial \Gamma'''}{\partial \mathbf c}\]where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.
- dkappa_by_dcoeff_vjp(v)
This function returns the vector Jacobian product
\[v^T \frac{\partial \kappa}{\partial \mathbf{c}}\]where \(\mathbf{c}\) are the curve dofs and \(\kappa\) is the curvature.
- dtorsion_by_dcoeff_vjp(v)
This function returns the vector Jacobian product
\[v^T \frac{\partial \tau}{\partial \mathbf{c}}\]where \(\mathbf{c}\) are the curve dofs, and \(\tau\) is the torsion.
- gamma_impl(gamma, quadpoints)
This function returns the x,y,z coordinates of the curve \(\Gamma\).
- gammadash_impl(gammadash)
This function returns \(\Gamma'(\varphi)\), where \(\Gamma\) are the x, y, z coordinates of the curve.
- gammadashdash_impl(gammadashdash)
This function returns \(\Gamma''(\varphi)\), and \(\Gamma\) are the x, y, z coordinates of the curve.
- gammadashdashdash_impl(gammadashdashdash)
This function returns \(\Gamma'''(\varphi)\), and \(\Gamma\) are the x, y, z coordinates of the curve.
- set_dofs(self: simsoptpp.Curve, arg0: List[float]) None
- class simsopt.geo.curve.RotatedCurve(curve, phi, flip)
Bases:
simsoptpp.Curve
,simsopt.geo.curve.Curve
RotatedCurve inherits from the Curve base class. It takes an input a Curve, rotates it about the
z
axis by a toroidal anglephi
, and optionally completes a reflection whenflip=True
.- dgamma_by_dcoeff_impl(dgamma_by_dcoeff)
This function returns
\[\frac{\partial \Gamma}{\partial \mathbf c}\]where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.
- dgamma_by_dcoeff_vjp(v)
This function returns the vector Jacobian product
\[v^T \frac{\partial \Gamma}{\partial \mathbf c}\]where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.
- dgammadash_by_dcoeff_impl(dgammadash_by_dcoeff)
This function returns
\[\frac{\partial \Gamma'}{\partial \mathbf c}\]where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.
- dgammadash_by_dcoeff_vjp(v)
This function returns the vector Jacobian product
\[v^T \frac{\partial \Gamma'}{\partial \mathbf c}\]where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.
- dgammadashdash_by_dcoeff_impl(dgammadashdash_by_dcoeff)
This function returns
\[\frac{\partial \Gamma''}{\partial \mathbf c}\]where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.
- dgammadashdash_by_dcoeff_vjp(v)
This function returns the vector Jacobian product
\[v^T \frac{\partial \Gamma''}{\partial \mathbf c}\]where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.
- dgammadashdashdash_by_dcoeff_impl(dgammadashdashdash_by_dcoeff)
This function returns
\[\frac{\partial \Gamma'''}{\partial \mathbf c}\]where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.
- dgammadashdashdash_by_dcoeff_vjp(v)
This function returns the vector Jacobian product
\[v^T \frac{\partial \Gamma'''}{\partial \mathbf c}\]where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.
- gamma_impl(gamma, quadpoints)
This function returns the x,y,z coordinates of the curve, \(\Gamma\), where \(\Gamma\) are the x, y, z coordinates of the curve.
- gammadash_impl(gammadash)
This function returns \(\Gamma'(\varphi)\), where \(\Gamma\) are the x, y, z coordinates of the curve.
- gammadashdash_impl(gammadashdash)
This function returns \(\Gamma''(\varphi)\), where \(\Gamma\) are the x, y, z coordinates of the curve.
- gammadashdashdash_impl(gammadashdashdash)
This function returns \(\Gamma'''(\varphi)\), where \(\Gamma\) are the x, y, z coordinates of the curve.
- get_dofs()
RotatedCurve does not have any dofs of its own. This function returns null array
- num_dofs()
This function returns the number of dofs associated to the curve.
- set_dofs_impl(d)
RotatedCurve does not have any dofs of its own. This function does nothing.
- simsopt.geo.curve.create_equally_spaced_curves(ncurves, nfp, stellsym, R0=1.0, R1=0.5, order=6, numquadpoints=None)
Create
ncurves
curves of typeCurveXYZFourier
of orderorder
that will result in circular equally spaced coils (major radiusR0
and minor radiusR1
) after applyingcoils_via_symmetries
.Usage example: create 4 base curves, which are then rotated 3 times and flipped for stellarator symmetry:
base_curves = create_equally_spaced_curves(4, 3, stellsym=True) base_currents = [Current(1e5) for c in base_curves] coils = coils_via_symmetries(base_curves, base_currents, 3, stellsym=True)
- simsopt.geo.curve.curves_to_vtk(curves, filename)
- simsopt.geo.curve.incremental_arclength_pure(d1gamma)
This function is used in a Python+Jax implementation of the curve arc length formula.
- simsopt.geo.curve.incremental_arclength_vjp(d1gamma, v)
- simsopt.geo.curve.kappa_pure(d1gamma, d2gamma)
This function is used in a Python+Jax implementation of formula for curvature.
- simsopt.geo.curve.kappagrad0(d1gamma, d2gamma)
- simsopt.geo.curve.kappagrad1(d1gamma, d2gamma)
- simsopt.geo.curve.kappavjp0(d1gamma, d2gamma, v)
- simsopt.geo.curve.kappavjp1(d1gamma, d2gamma, v)
- simsopt.geo.curve.torsion_pure(d1gamma, d2gamma, d3gamma)
This function is used in a Python+Jax implementation of formula for torsion.
- simsopt.geo.curve.torsionvjp0(d1gamma, d2gamma, d3gamma, v)
- simsopt.geo.curve.torsionvjp1(d1gamma, d2gamma, d3gamma, v)
- simsopt.geo.curve.torsionvjp2(d1gamma, d2gamma, d3gamma, v)
simsopt.geo.curvehelical module
- class simsopt.geo.curvehelical.CurveHelical(quadpoints, order, n0, l0, R0, r0)
Bases:
simsopt.geo.curve.JaxCurve
Curve representation of a helical coil. The helical coil positions are specified by a poloidal angle eta that is a function of the toroidal angle phi with Fourier coefficients \(A_k\) and \(B_k\). The poloidal angle is represented as
\[\eta = m_0 \phi/l_0 + \sum_k A_k \cos(n_0 \phi k/l_0) + B_k \sin(n_0 \phi k/l_0)\]- Parameters
quadpoints – number of grid points/resolution along the curve;
order – number of fourier coefficients, i.e., the length of the array A (or B);
n0 – toroidal periodicity/number of field periods
l0 – number of \(2\pi\) turns in \(\phi\)
R0 – major radius
r0 – minor radius
- get_dofs(self: simsoptpp.Curve) List[float]
- num_dofs(self: simsoptpp.Curve) int
- set_dofs_impl(self: simsoptpp.Curve, arg0: List[float]) None
- simsopt.geo.curvehelical.jaxHelicalfouriercurve_pure(dofs, quadpoints, order, n0, l0, R0, r0)
simsopt.geo.curveobjectives module
- class simsopt.geo.curveobjectives.CurveLength(curve)
Bases:
simsopt._core.graph_optimizable.Optimizable
CurveLength is a class that computes the length of a curve, i.e.
\[J = \int_{\text{curve}}~dl.\]- J()
This returns the value of the quantity.
- dJ(*args, partials=False, **kwargs)
- return_fn_map: Dict[str, Callable] = {'J': <function CurveLength.J>, 'dJ': <function derivative_dec.<locals>._derivative_dec>}
- class simsopt.geo.curveobjectives.LpCurveCurvature(curve, p, desired_length=None)
Bases:
simsopt._core.graph_optimizable.Optimizable
This class computes a penalty term based on the \(L_p\) norm of the curve’s curvature, and penalizes where the local curve curvature exceeds a threshold
\[J = \frac{1}{p} \int_{\text{curve}} \text{max}(\kappa - \kappa_0, 0)^p ~dl\]where \(\kappa_0\) is a threshold curvature. When
desired_length
is provided, \(\kappa_0 = 2 \pi / \text{desired_length}\), otherwise \(\kappa_0=0\).- J()
This returns the value of the quantity.
- dJ(*args, partials=False, **kwargs)
- return_fn_map: Dict[str, Callable] = {'J': <function LpCurveCurvature.J>, 'dJ': <function derivative_dec.<locals>._derivative_dec>}
- class simsopt.geo.curveobjectives.LpCurveTorsion(curve, p)
Bases:
simsopt._core.graph_optimizable.Optimizable
LpCurveTorsion is a class that computes a penalty term based on the \(L_p\) norm of the curve’s torsion:
\[J = \frac{1}{p} \int_{\text{curve}} \tau^p ~dl.\]- J()
This returns the value of the quantity.
- dJ(*args, partials=False, **kwargs)
- return_fn_map: Dict[str, Callable] = {'J': <function LpCurveTorsion.J>, 'dJ': <function derivative_dec.<locals>._derivative_dec>}
- simsopt.geo.curveobjectives.Lp_curvature_pure(kappa, gammadash, p, desired_kappa)
This function is used in a Python+Jax implementation of the curvature penalty term.
- simsopt.geo.curveobjectives.Lp_torsion_pure(torsion, gammadash, p)
This function is used in a Python+Jax implementation of the formula for the torsion penalty term.
- class simsopt.geo.curveobjectives.MinimumDistance(curves, minimum_distance)
Bases:
simsopt._core.graph_optimizable.Optimizable
MinimumDistance is a class that computes
\[J = \sum_{i = 1}^{\text{num_coils}} \sum_{j = 1}^{i-1} d_{i,j}\]where
\[\begin{split}d_{i,j} = \int_{\text{curve}_i} \int_{\text{curve}_j} \max(0, d_{\min} - \| \mathbf{r}_i - \mathbf{r}_j \|_2)^2 ~dl_j ~dl_i\\\end{split}\]and \(\mathbf{r}_i\), \(\mathbf{r}_j\) are points on coils \(i\) and \(j\), respectively. \(d_\min\) is a desired threshold minimum intercoil distance. This penalty term is zero when the points on coil \(i\) and coil \(j\) lie more than \(d_\min\) away from one another, for \(i, j \in \{1, \cdots, \text{num_coils}\}\)
- J()
This returns the value of the quantity.
- dJ(*args, partials=False, **kwargs)
- 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 MinimumDistance.J>, 'dJ': <function derivative_dec.<locals>._derivative_dec>}
- simsopt.geo.curveobjectives.curve_length_pure(l)
This function is used in a Python+Jax implementation of the curve length formula.
- simsopt.geo.curveobjectives.distance_pure(gamma1, l1, gamma2, l2, minimum_distance)
This function is used in a Python+Jax implementation of the distance formula.
simsopt.geo.curverzfourier module
- class simsopt.geo.curverzfourier.CurveRZFourier(quadpoints, order, nfp, stellsym)
Bases:
simsoptpp.CurveRZFourier
,simsopt.geo.curve.Curve
- CurveRZFourier is a curve that is represented in cylindrical
coordinates using the following Fourier series:
\[\begin{split}r(\phi) &= \sum_{m=0}^{\text{order}} x_{c,m}\cos(n_{\text{fp}} m \phi) + \sum_{m=1}^{\text{order}} x_{s,m}\sin(n_{\text{fp}} m \phi) \\ z(\phi) &= \sum_{m=0}^{\text{order}} z_{c,m}\cos(n_{\text{fp}} m \phi) + \sum_{m=1}^{\text{order}} z_{s,m}\sin(n_{\text{fp}} m \phi)\end{split}\]If
stellsym = True
, then the \(\sin\) terms for \(r\) and the \(\cos\) terms for \(z\) are zero.For the
stellsym = False
case, the dofs are stored in the order\[[r_{c,0}, \cdots, r_{c,\text{order}}, r_{s,1}, \cdots, r_{s,\text{order}}, z_{c,0},....]\]or in the
stellsym = True
case they are stored\[[r_{c,0},...,r_{c,order},z_{s,1},...,z_{s,order}]\]
- get_dofs()
This function returns the dofs associated to this object.
- set_dofs(dofs)
This function sets the dofs associated to this object.
simsopt.geo.curvexyzfourier module
- class simsopt.geo.curvexyzfourier.CurveXYZFourier(quadpoints, order)
Bases:
simsoptpp.CurveXYZFourier
,simsopt.geo.curve.Curve
CurveXYZFourier is a curve that is represented in Cartesian coordinates using the following Fourier series:
\[\begin{split}x(\phi) &= \sum_{m=0}^{\text{order}} x_{c,m}\cos(m\phi) + \sum_{m=1}^{\text{order}} x_{s,m}\sin(m\phi) \\ y(\phi) &= \sum_{m=0}^{\text{order}} y_{c,m}\cos(m\phi) + \sum_{m=1}^{\text{order}} y_{s,m}\sin(m\phi) \\ z(\phi) &= \sum_{m=0}^{\text{order}} z_{c,m}\cos(m\phi) + \sum_{m=1}^{\text{order}} z_{s,m}\sin(m\phi)\end{split}\]The dofs are stored in the order
\[[x_{c,0}, x_{s,1}, x_{c,1},\cdots x_{s,\text{order}}, x_{c,\text{order}},y_{c,0},y_{s,1},y_{c,1},\cdots]\]- get_dofs()
This function returns the dofs associated to this object.
- static load_curves_from_file(filename, order=None, ppp=20, delimiter=',')
This function loads a file containing Fourier coefficients for several coils. The file is expected to have
6*num_coils
many columns, andorder+1
many rows. The columns are in the following order,sin_x_coil1, cos_x_coil1, sin_y_coil1, cos_y_coil1, sin_z_coil1, cos_z_coil1, sin_x_coil2, cos_x_coil2, sin_y_coil2, cos_y_coil2, sin_z_coil2, cos_z_coil2, …
- set_dofs(dofs)
This function sets the dofs associated to this object.
- class simsopt.geo.curvexyzfourier.JaxCurveXYZFourier(quadpoints, order)
Bases:
simsopt.geo.curve.JaxCurve
A Python+Jax implementation of the CurveXYZFourier class. There is actually no reason why one should use this over the C++ implementation in
simsoptpp
, but the point of this class is to illustrate how jax can be used to define a geometric object class and calculate all the derivatives (both with respect to dofs and with respect to the angle \(\phi\)) automatically.- get_dofs()
This function returns the dofs associated to this object.
- num_dofs()
This function returns the number of dofs associated to this object.
- set_dofs_impl(dofs)
This function sets the dofs associated to this object.
- simsopt.geo.curvexyzfourier.jaxfouriercurve_pure(dofs, quadpoints, order)
simsopt.geo.jit module
- simsopt.geo.jit.jit(fun)
simsopt.geo.plot module
- simsopt.geo.plot.fix_matplotlib_3d(ax)
Make axes of 3D plot have equal scale so that spheres appear as spheres, cubes as cubes, etc.. This is one possible solution to Matplotlib’s
ax.set_aspect('equal')
andax.axis('equal')
not working for 3D. This function is to be called after objects have been plotted.This function was taken from https://stackoverflow.com/questions/13685386/matplotlib-equal-unit-length-with-equal-aspect-ratio-z-axis-is-not-equal-to
- Parameters
ax – a matplotlib axis, e.g., as output from
plt.gca()
.
simsopt.geo.qfmsurface module
- class simsopt.geo.qfmsurface.QfmSurface(biotsavart, surface, label, targetlabel)
Bases:
object
QfmSurface is used to compute a quadratic-flux minimizing surface, defined as the minimum of the objective function,
\[f(S) = \frac{\int_{S} d^2 x \, \left(\textbf{B} \cdot \hat{\textbf{n}}\right)^2}{\int_{S} d^2 x \, B^2}\]subject to a constraint on the surface label, such as the area, volume, or toroidal flux. This objective, \(f\), is computed in
QfmResidual
(seesurfaceobjectives.py
). If magnetic surfaces exists, \(f\) will be zero. If not, QFM surfaces approximate flux surfaces.The label constraint can be enforced with a penalty formulation, defined in
qfm_penalty_constraints()
, and whose minimium is computed with LBFGS-B by :func:`minimize_qfm_penalty_constraints_LBFGS()’.Alternatively, the label constraint can be enforced with a constrained optimization algorithm. This constrained optimization problem is solved with the SLSQP algorithm by
minimize_qfm_exact_constraints_SLSQP()
.- minimize_qfm(tol=0.001, maxiter=1000, method='SLSQP', constraint_weight=1.0)
This function tries to find the surface that approximately solves
\[\min_{S} f(S)\]subject to
\[\texttt{label} = \texttt{labeltarget}\]where \(f(S)\) is the QFM residual. This is done using method, either ‘SLSQP’ or ‘LBFGS’. If ‘LBFGS’ is chosen, a penalized formulation is used defined by constraint_weight. If ‘SLSQP’ is chosen, constraint_weight is not used, and the constrained optimization problem is solved.
- minimize_qfm_exact_constraints_SLSQP(tol=0.001, maxiter=1000)
This function tries to find the surface that approximately solves
\[\min_{S} f(S)\]subject to
\[\texttt{label} = \texttt{labeltarget}\]where \(f(S)\) is the QFM residual. This is done using SLSQP.
- minimize_qfm_penalty_constraints_LBFGS(tol=0.001, maxiter=1000, constraint_weight=1.0)
This function tries to find the surface that approximately solves
\[\min_S f(S) + 0.5 \texttt{constraint_weight} (\texttt{label} - \texttt{labeltarget})^2\]where \(f(S)\) is the QFM residual defined in
surfaceobjectives.py.
This is done using LBFGS.
- qfm_label_constraint(x, derivatives=0)
Returns the residual
\[0.5 (\texttt{label} - \texttt{labeltarget})^2\]and optionally the gradient wrt surface params.
- qfm_objective(x, derivatives=0)
Returns the residual
\[f(S)\]and optionally the gradient wrt surface params where \(f(S)\) is the QFM residual defined in surfaceobjectives.py.
- qfm_penalty_constraints(x, derivatives=0, constraint_weight=1)
Returns the residual
\[f(S) + 0.5 \texttt{constraint_weight} (\texttt{label} - \texttt{labeltarget})^2\]and optionally the gradient and Hessian wrt surface params where \(f(S)\) is the QFM residual defined in
surfaceobjectives.py
.
simsopt.geo.surface module
- class simsopt.geo.surface.Surface(**kwargs)
Bases:
simsopt._core.graph_optimizable.Optimizable
Surface
is a base class for various representations of toroidal surfaces in simsopt.A
Surface
is modelled as a function \(\Gamma:[0, 1] \times [0, 1] \to R^3\) and is evaluated at quadrature points \(\{\phi_1, \ldots, \phi_{n_\phi}\}\times\{\theta_1, \ldots, \theta_{n_\theta}\}\).- RANGE_FIELD_PERIOD = 'field period'
- RANGE_FULL_TORUS = 'full torus'
- RANGE_HALF_PERIOD = 'half period'
- arclength_poloidal_angle()
Computes poloidal angle based on arclenth along magnetic surface at constant phi. The resulting angle is in the range [0,1]. This is required for evaluating the adjoint shape gradient for free-boundary calculations.
- Returns
- 2d array (numquadpoints_phi,numquadpoints_theta)
of arclength poloidal angle
- Return type
theta_arclength
- aspect_ratio()
Note: cylindrical coordinates are \((R, \phi, Z)\), where \(\phi \in [-\pi,\pi)\) and the angles that parametrize the surface are \((\varphi, \theta) \in [0,1)^2\) For a given surface, this function computes its aspect ratio using the VMEC definition:
\[AR = R_{\text{major}} / R_{\text{minor}}\]where
\[\begin{split}R_{\text{minor}} &= \sqrt{ \overline{A} / \pi } \\ R_{\text{major}} &= \frac{V}{2 \pi^2 R_{\text{minor}}^2}\end{split}\]and \(V\) is the volume enclosed by the surface, and \(\overline{A}\) is the average cross sectional area. The main difficult part of this calculation is the mean cross sectional area. This is given by the integral
\[\overline{A} = \frac{1}{2\pi} \int_{S_{\phi}} ~dS ~d\phi\]where \(S_\phi\) is the cross section of the surface at the cylindrical angle \(\phi\). Note that \(\int_{S_\phi} ~dS\) can be rewritten as a line integral
\[\begin{split}\int_{S_\phi}~dS &= \int_{S_\phi} ~dR dZ \\ &= \int_{\partial S_\phi} [R,0] \cdot \mathbf n/\|\mathbf n\| ~dl \\ &= \int^1_{0} R \frac{\partial Z}{\partial \theta}~d\theta\end{split}\]where \(\mathbf n = [n_R, n_Z] = [\partial Z/\partial \theta, -\partial R/\partial \theta]\) is the outward pointing normal.
Consider the surface in cylindrical coordinates terms of its angles \([R(\varphi,\theta), \phi(\varphi,\theta), Z(\varphi,\theta)]\). The boundary of the cross section \(\partial S_\phi\) is given by the points \(\theta\rightarrow[R(\varphi(\phi,\theta),\theta),\phi, Z(\varphi(\phi,\theta),\theta)]\) for fixed \(\phi\). The cross sectional area of \(S_\phi\) becomes
\[\int^{1}_{0} R(\varphi(\phi,\theta),\theta) \frac{\partial}{\partial \theta}[Z(\varphi(\phi,\theta),\theta)] ~d\theta\]Now, substituting this into the formula for the mean cross sectional area, we have
\[\overline{A} = \frac{1}{2\pi}\int^{\pi}_{-\pi}\int^{1}_{0} R(\varphi(\phi,\theta),\theta) \frac{\partial}{\partial \theta}[Z(\varphi(\phi,\theta),\theta)] ~d\theta ~d\phi\]Instead of integrating over cylindrical \(\phi\), let’s complete the change of variables and integrate over \(\varphi\) using the mapping:
\[[\phi,\theta] \leftarrow [\text{atan2}(y(\varphi,\theta), x(\varphi,\theta)), \theta]\]After the change of variables, the integral becomes:
\[\overline{A} = \frac{1}{2\pi}\int^{1}_{0}\int^{1}_{0} R(\varphi,\theta) \left[\frac{\partial Z}{\partial \varphi} \frac{\partial \varphi}{d \theta} + \frac{\partial Z}{\partial \theta} \right] \text{det} J ~d\theta ~d\varphi\]where \(\text{det}J\) is the determinant of the mapping’s Jacobian.
- cross_section(phi, thetas=None)
This function takes in a cylindrical angle \(\phi\) and returns the cross section of the surface in that plane evaluated at thetas. This is done using the method of bisection. This function takes in a cylindrical angle \(\phi\) and returns the cross section of the surface in that plane evaluated at thetas. This is done using the method of bisection.
This function assumes that the surface intersection with the plane is a single curve.
- get_quadpoints(quadpoints_theta=None, range='full torus', nphi=None, ntheta=None, nfp=1)
This function is used to set the theta and phi grid points for Surface subclasses. It is typically called in the constructor of each Surface subclass.
For more information about the arguments
nphi
,ntheta
,range
,quadpoints_phi
, andquadpoints_theta
, see the general documentation on Surfaces.- Parameters
nfp – The number of field periods.
nphi – Number of grid points \(\phi_j\) in the toroidal angle \(\phi\).
ntheta – Number of grid points \(\theta_j\) in the toroidal angle \(\theta\).
range – Toroidal extent of the \(\phi\) grid. Set to
"full torus"
(or equivalentlySurface.RANGE_FULL_TORUS
) to generate points up to 1 (with no point at 1). Set to"field period"
(or equivalentlySurface.RANGE_FIELD_PERIOD
) to generate points up to \(1/n_{fp}\) (with no point at \(1/n_{fp}\)). Set to"half period"
(or equivalentlySurface.RANGE_HALF_PERIOD
) to generate points up to \(1/(2 n_{fp})\) (with no point at \(1/(2 n_{fp})\)). Ifquadpoints_phi
is specified,range
is irrelevant.quadpoints_phi – Set this to a list or 1D array to set the \(\phi_j\) grid points directly.
quadpoints_theta – Set this to a list or 1D array to set the \(\theta_j\) grid points directly.
- interpolate_on_arclength_grid(function, theta_evaluate)
Interpolate function onto the theta_evaluate grid in the arclength poloidal angle. This is required for evaluating the adjoint shape gradient for free-boundary calculations.
- Returns
- 2d array (numquadpoints_phi,numquadpoints_theta)
defining interpolated function on arclength angle along curve at constant phi
- Return type
function_interpolated
- plot(engine='matplotlib', ax=None, show=True, close=False, axis_equal=True, plot_normal=False, plot_derivative=False, wireframe=True, **kwargs)
Plot the surface in 3D using matplotlib/mayavi/plotly.
- Parameters
engine – Selects the graphics engine. Currently supported options are
"matplotlib"
(default),"mayavi"
, and"plotly"
.ax – The figure/axis to be plotted on. This argument is useful when plotting multiple objects on the same axes. If equal to the default
None
, a new axis will be created.show – Whether to call the
show()
function of the graphics engine. Should be set toFalse
if more objects will be plotted on the same axes.close – Whether to close the seams in the surface where the angles jump back to 0.
axis_equal – For matplotlib, whether to adjust the scales of the x, y, and z axes so distances in each direction appear equal.
plot_normal – Whether to plot the surface normal vectors. Only implemented for mayavi.
plot_derivative – Whether to plot the surface derivatives. Only implemented for mayavi.
wireframe – Whether to plot the wireframe in Mayavi.
kwargs – Any additional arguments to pass to the plotting function, like
color='r'
.
Note: the
ax
andshow
parameters can be used to plot more than one surface:ax = surface1.plot(show=False) ax = surface2.plot(ax=ax, show=False) surface3.plot(ax=ax, show=True)
- Returns
An axis which could be passed to a further call to the graphics engine so multiple objects are shown together.
- abstract to_RZFourier()
Return a
simsopt.geo.surfacerzfourier.SurfaceRZFourier
instance corresponding to the shape of this surface. All subclasses should implement this abstract method.
- to_vtk(filename, extra_data=None)
- class simsopt.geo.surface.SurfaceClassifier(surface, p=1, h=0.05)
Bases:
object
Takes in a toroidal surface and constructs an interpolant of the signed distance function \(f:R^3\to R\) that is positive inside the volume contained by the surface, (approximately) zero on the surface, and negative outisde the volume contained by the surface.
- __init__(surface, p=1, h=0.05)
- Parameters
surface – the surface to contruct the distance from.
p – degree of the interpolant
h – grid resolution of the interpolant
- evaluate(xyz)
- to_vtk(filename, h=0.01)
- simsopt.geo.surface.signed_distance_from_surface(xyz, surface)
Compute the signed distances from points
xyz
to a surface. The sign is positive for points inside the volume surrounded by the surface.
simsopt.geo.surfacegarabedian module
- class simsopt.geo.surfacegarabedian.SurfaceGarabedian(nfp=1, mmax=1, mmin=0, nmax=0, nmin=None, nphi=None, ntheta=None, range='full torus', quadpoints_phi=None, quadpoints_theta=None)
Bases:
simsoptpp.Surface
,simsopt.geo.surface.Surface
SurfaceGarabedian
represents a toroidal surface for which the shape is parameterized using Garabedian’s \(\Delta_{m,n}\) coefficients:\[R + i Z = e^{i u} \sum_{m = m_\min}^{m_\max} \sum_{n = n_\min}^{n_\max} \Delta_{m,n} e^{-i m u + i n v}\]where \(u = 2 \pi \theta\) is a poloidal angle on \([0, 2\pi]\), and \(v\) is the standard toroidal angle on \([0, 2\pi]\).
The present implementation assumes stellarator symmetry. Note that non-stellarator-symmetric surfaces require that the \(\Delta_{m,n}\) coefficients be imaginary.
For more information about the arguments
nphi
,ntheta
,range
,quadpoints_phi
, andquadpoints_theta
, see the general documentation on Surfaces.- Parameters
nfp – The number of field periods.
mmin – Minimum poloidal mode number \(m\) included (usually 0 or negative).
mmax – Maximum poloidal mode number \(m\) included.
nmin – Minimum toroidal mode number \(n\) included (usually negative). If
None
,nmin = -nmax
will be used.nmax – Maximum toroidal mode number \(n\) included.
nphi – Number of grid points \(\phi_j\) in the toroidal angle \(\phi\).
ntheta – Number of grid points \(\theta_j\) in the toroidal angle \(\theta\).
range – Toroidal extent of the \(\phi\) grid. Set to
"full torus"
(or equivalentlySurfaceGarabedian.RANGE_FULL_TORUS
) to generate points up to 1 (with no point at 1). Set to"field period"
(or equivalentlySurfaceGarabedian.RANGE_FIELD_PERIOD
) to generate points up to \(1/n_{fp}\) (with no point at \(1/n_{fp}\)). Set to"half period"
(or equivalentlySurfaceGarabedian.RANGE_HALF_PERIOD
) to generate points up to \(1/(2 n_{fp})\) (with no point at \(1/(2 n_{fp})\)). Ifquadpoints_phi
is specified,range
is irrelevant.quadpoints_phi – Set this to a list or 1D array to set the \(\phi_j\) grid points directly.
quadpoints_theta – Set this to a list or 1D array to set the \(\theta_j\) grid points directly.
- property Delta
- area()
Return the area of the surface.
- area_volume()
Compute the surface area and the volume enclosed by the surface.
- fix_range(mmin, mmax, nmin, nmax, fixed=True)
Fix the DOFs for a range of m and n values.
All modes with m in the interval [mmin, mmax] and n in the interval [nmin, nmax] will have their fixed property set to the value of the ‘fixed’ parameter. Note that mmax and nmax are included (unlike the upper bound in python’s range(min, max).)
- classmethod from_RZFourier(surf)
Create a SurfaceGarabedian from a SurfaceRZFourier object of the identical shape.
For a derivation of the transformation here, see https://terpconnect.umd.edu/~mattland/assets/notes/toroidal_surface_parameterizations.pdf
- get_Delta(m, n)
Return a particular \(\Delta_{m,n}\) coefficient.
- get_dofs()
Return a 1D numpy array with all the degrees of freedom.
- return_fn_map: Dict[str, Callable] = {'area': <function SurfaceGarabedian.area>, 'aspect-ratio': <function Surface.aspect_ratio>, 'volume': <function SurfaceGarabedian.volume>}
- set_Delta(m, n, val)
Set a particular \(\Delta_{m,n}\) coefficient.
- set_dofs(x)
Set the shape coefficients from a 1D list/array
- to_RZFourier()
Return a SurfaceRZFourier object with the identical shape.
For a derivation of the transformation here, see https://terpconnect.umd.edu/~mattland/assets/notes/toroidal_surface_parameterizations.pdf
- volume()
Return the volume of the surface.
simsopt.geo.surfacehenneberg module
- class simsopt.geo.surfacehenneberg.SurfaceHenneberg(nfp: int = 1, alpha_fac: int = 1, mmax: int = 1, nmax: int = 0, nphi: Optional[int] = None, ntheta: Optional[int] = None, range: str = 'full torus', quadpoints_phi: Optional[Union[Sequence[numbers.Real], nptyping.types._ndarray.NDArray[None, nptyping.types._number.Float[float, numpy.floating]]]] = None, quadpoints_theta: Optional[Union[Sequence[numbers.Real], nptyping.types._ndarray.NDArray[None, nptyping.types._number.Float[float, numpy.floating]]]] = None)
Bases:
simsoptpp.Surface
,simsopt.geo.surface.Surface
This class represents a toroidal surface using the parameterization in Henneberg, Helander, and Drevlak, Journal of Plasma Physics 87, 905870503 (2021). The main benefit of this representation is that there is no freedom in the poloidal angle, i.e. \(\theta\) is uniquely defined, in contrast to other parameterizations like
SurfaceRZFourier
. Stellarator symmetry is assumed.In this representation by Henneberg et al, the cylindrical coordinates \((R,\phi,Z)\) are written in terms of a unique poloidal angle \(\theta\) as follows:
\[\begin{split}R(\theta,\phi) = R_0^H(\phi) + \rho(\theta,\phi) \cos(\alpha\phi) - \zeta(\theta,\phi) \sin(\alpha\phi), \\ Z(\theta,\phi) = Z_0^H(\phi) + \rho(\theta,\phi) \sin(\alpha\phi) + \zeta(\theta,\phi) \cos(\alpha\phi),\end{split}\]where
\[\begin{split}R_0^H(\phi) &=& \sum_{n=0}^{nmax} R_{0,n}^H \cos(n_{fp} n \phi), \\ Z_0^H(\phi) &=& \sum_{n=1}^{nmax} Z_{0,n}^H \sin(n_{fp} n \phi), \\ \zeta(\theta,\phi) &=& \sum_{n=0}^{nmax} b_n \cos(n_{fp} n \phi) \sin(\theta - \alpha \phi), \\ \rho(\theta,\phi) &=& \sum_{n,m} \rho_{n,m} \cos(m \theta - n_{fp} n \phi - \alpha \phi).\end{split}\]The continuous degrees of freedom are \(\{\rho_{m,n}, b_n, R_{0,n}^H, Z_{0,n}^H\}\). These variables correspond to the attributes
rhomn
,bn
,R0nH
, andZ0nH
respectively, which are all numpy arrays. There is also a discrete degree of freedom \(\alpha\) which should be \(\pm n_{fp}/2\) where \(n_{fp}\) is the number of field periods. The attributealpha_fac
corresponds to \(2\alpha/n_{fp}\), soalpha_fac
is either 1, 0, or -1. Usingalpha_fac = 0
is appropriate for axisymmetry, while values of 1 or -1 are appropriate for a stellarator, depending on the handedness of the rotating elongation.For \(R_{0,n}^H\) and \(b_n\), \(n\) is 0 or any positive integer up through
nmax
(inclusive). For \(Z_{0,n}^H\), \(n\) is any positive integer up throughnmax
. For \(\rho_{m,n}\), \(m\) is an integer from 0 throughmmax
(inclusive). For positive values of \(m\), \(n\) can be any integer from-nmax
throughnmax
. For \(m=0\), \(n\) is restricted to integers from 1 throughnmax
. Note that we exclude the element of \(\rho_{m,n}\) with \(m=n=0\), because this degree of freedom is already represented in \(R_{0,0}^H\).For the 2D array
rhomn
, functionsset_rhomn()
andget_rhomn()
are provided for convenience so you can specifyn
, since the corresponding array index is shifted bynmax
. There are no corresponding functions for the 1D arraysR0nH
,Z0nH
, andbn
since these arrays all have a first index corresponding ton=0
.For more information about the arguments
nphi
,ntheta
,range
,quadpoints_phi
, andquadpoints_theta
, see the general documentation on Surfaces.- Parameters
nfp – The number of field periods.
alpha_fac – Should be +1 or -1 for a stellarator, depending on the handedness by which the elongation rotates, or 0 for axisymmetry.
mmax – Maximum poloidal mode number included.
nmax – Maximum toroidal mode number included, divided by
nfp
.nphi – Number of grid points \(\phi_j\) in the toroidal angle \(\phi\).
ntheta – Number of grid points \(\theta_j\) in the toroidal angle \(\theta\).
range – Toroidal extent of the \(\phi\) grid. Set to
"full torus"
(or equivalentlySurfaceHenneberg.RANGE_FULL_TORUS
) to generate points up to 1 (with no point at 1). Set to"field period"
(or equivalentlySurfaceHenneberg.RANGE_FIELD_PERIOD
) to generate points up to \(1/n_{fp}\) (with no point at \(1/n_{fp}\)). Set to"half period"
(or equivalentlySurfaceHenneberg.RANGE_HALF_PERIOD
) to generate points up to \(1/(2 n_{fp})\) (with no point at \(1/(2 n_{fp})\)). Ifquadpoints_phi
is specified,range
is irrelevant.quadpoints_phi – Set this to a list or 1D array to set the \(\phi_j\) grid points directly.
quadpoints_theta – Set this to a list or 1D array to set the \(\theta_j\) grid points directly.
- _validate_mn(m, n)
Check whether given (m, n) values are allowed for \(\rho_{m,n}\).
- allocate()
Create the arrays for the continuous degrees of freedom. Also set the names of the dofs.
- fixed_range(mmax, nmax, fixed=True)
Set the
fixed
property for a range ofm
andn
values.All modes with
m <= mmax
and|n| <= nmax
will have their fixed property set to the value of thefixed
parameter. Note thatmmax
andnmax
are included.Both
mmax
andnmax
must be >= 0.For any value of
mmax
, thefixed
properties ofR0nH
,Z0nH
, andrhomn
are set. Thefixed
properties ofbn
are set only ifmmax > 0
. In other words, thebn
modes are treated as havingm=1
.
- classmethod from_RZFourier(surf, alpha_fac: int, mmax: Optional[int] = None, nmax: Optional[int] = None, ntheta: Optional[int] = None, nphi: Optional[int] = None)
Convert a
SurfaceRZFourier
surface to aSurfaceHenneberg
surface.- Parameters
surf – The
SurfaceRZFourier
object to convert.mmax – Maximum poloidal mode number to include in the new surface. If
None
, the valuempol
from the old surface will be used.nmax – Maximum toroidal mode number to include in the new surface. If
None
, the valuentor
from the old surface will be used.ntheta – Number of grid points in the poloidal angle used for the transformation. If
None
, the value3 * ntheta
will be used.nphi – Number of grid points in the toroidal angle used for the transformation. If
None
, the value3 * nphi
will be used.
- gamma_impl(data, quadpoints_phi, quadpoints_theta)
Evaluate the position vector on the surface in Cartesian coordinates, for a tensor product grid of points in theta and phi.
- gamma_lin(data, quadpoints_phi, quadpoints_theta)
Evaluate the position vector on the surface in Cartesian coordinates, for a list of (phi, theta) points.
- gammadash1_impl(data)
Evaluate the derivative of the position vector with respect to the toroidal angle phi.
- gammadash2_impl(data)
Evaluate the derivative of the position vector with respect to theta.
- get_dofs()
Return a 1D numpy array with all the degrees of freedom.
- get_rhomn(m, n)
Return a particular \(\rho_{m,n}\) coefficient.
- num_dofs()
Return the number of degrees of freedom.
- set_dofs(self: simsoptpp.Surface, arg0: List[float]) None
- set_dofs_impl(v)
Set the shape coefficients from a 1D list/array
- set_rhomn(m, n, val)
Set a particular \(\rho_{m,n}\) coefficient.
- to_RZFourier()
Return a
SurfaceRZFourier
object with the identical shape. This routine implements eq (4.5)-(4.6) in the Henneberg paper, plus m=0 terms for R0 and Z0.
simsopt.geo.surfaceobjectives module
- class simsopt.geo.surfaceobjectives.Area(surface)
Bases:
simsopt._core.graph_optimizable.Optimizable
Wrapper class for surface area label.
- J()
Compute the area of a surface.
- d2J_by_dsurfacecoefficientsdsurfacecoefficients()
Calculate the second derivatives with respect to the surface coefficients.
- dJ_by_dsurfacecoefficients()
Calculate the derivatives with respect to the surface coefficients.
- class simsopt.geo.surfaceobjectives.QfmResidual(surface, biotsavart)
Bases:
simsopt._core.graph_optimizable.Optimizable
For a given surface \(S\), this class computes the residual
\[f(S) = \frac{\int_{S} d^2 x \, (\textbf{B} \cdot \hat{\textbf{n}})^2}{\int_{S} d^2 x \, B^2}\]where \(\textbf{B}\) is the magnetic field from
biotsavart
, \(\hat{\textbf{n}}\) is the unit normal on a given surface, and the integration is performed over the surface. Derivatives are computed wrt the surface dofs.- J()
- dJ_by_dsurfacecoefficients()
Calculate the derivatives with respect to the surface coefficients
- invalidate_cache()
- recompute_bell(parent=None)
Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.
This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.
Need to be implemented by classes that provide a dof_setter for external handling of DOFs.
- class simsopt.geo.surfaceobjectives.ToroidalFlux(surface, biotsavart, idx=0)
Bases:
simsopt._core.graph_optimizable.Optimizable
Given a surface and Biot Savart kernel, this objective calculates
\[\begin{split}J &= \int_{S_{\varphi}} \mathbf{B} \cdot \mathbf{n} ~ds, \\ &= \int_{S_{\varphi}} \text{curl} \mathbf{A} \cdot \mathbf{n} ~ds, \\ &= \int_{\partial S_{\varphi}} \mathbf{A} \cdot \mathbf{t}~dl,\end{split}\]where \(S_{\varphi}\) is a surface of constant \(\varphi\), and \(\mathbf A\) is the magnetic vector potential.
- J()
Compute the toroidal flux on the surface where \(\varphi = \texttt{quadpoints_varphi}[\texttt{idx}]\).
- d2J_by_dsurfacecoefficientsdsurfacecoefficients()
Calculate the second derivatives with respect to the surface coefficients.
- dJ_by_dsurfacecoefficients()
Calculate the derivatives with respect to the surface coefficients.
- invalidate_cache()
- recompute_bell(parent=None)
Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.
This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.
Need to be implemented by classes that provide a dof_setter for external handling of DOFs.
- class simsopt.geo.surfaceobjectives.Volume(surface)
Bases:
simsopt._core.graph_optimizable.Optimizable
Wrapper class for volume label.
- J()
Compute the volume enclosed by the surface.
- d2J_by_dsurfacecoefficientsdsurfacecoefficients()
Calculate the second derivatives with respect to the surface coefficients.
- dJ_by_dsurfacecoefficients()
Calculate the derivatives with respect to the surface coefficients.
- simsopt.geo.surfaceobjectives.boozer_surface_residual(surface, iota, G, biotsavart, derivatives=0)
For a given surface, this function computes the residual
\[G\mathbf B_\text{BS}(\mathbf x) - ||\mathbf B_\text{BS}(\mathbf x)||^2 (\mathbf x_\varphi + \iota \mathbf x_\theta)\]as well as the derivatives of this residual with respect to surface dofs, iota, and G. In the above, \(\mathbf x\) are points on the surface, \(\iota\) is the rotational transform on that surface, and \(\mathbf B_{\text{BS}}\) is the magnetic field computed using the Biot-Savart law.
\(G\) is known for exact boozer surfaces, so if
G=None
is passed, then that value is used instead.
- simsopt.geo.surfaceobjectives.parameter_derivatives(surface: simsopt.geo.surface.Surface, shape_gradient: nptyping.types._ndarray.NDArray[typing.Any, typing.Any, nptyping.types._number.Float[float, numpy.floating]]) nptyping.types._ndarray.NDArray[typing.Any, nptyping.types._number.Float[float, numpy.floating]]
Converts the shape gradient of a given figure of merit, \(f\), to derivatives with respect to parameters defining a surface. For a perturbation to the surface \(\delta \vec{x}\), the resulting perturbation to the objective function is
\[\delta f(\delta \vec{x}) = \int d^2 x \, G \delta \vec{x} \cdot \vec{n}\]where \(G\) is the shape gradient and \(\vec{n}\) is the unit normal. Given \(G\), the parameter derivatives are then computed as
\[\frac{\partial f}{\partial \Omega} = \int d^2 x \, G \frac{\partial\vec{x}}{\partial \Omega} \cdot \vec{n},\]where \(\Omega\) is any parameter of the surface.
- Parameters
surface – The surface to use for the computation
shape_gradient – 2d array of size (numquadpoints_phi,numquadpoints_theta)
- Returns
1d array of size (ndofs)
simsopt.geo.surfacerzfourier module
- class simsopt.geo.surfacerzfourier.SurfaceRZFourier(nfp=1, stellsym=True, mpol=1, ntor=0, nphi=None, ntheta=None, range='full torus', quadpoints_phi=None, quadpoints_theta=None)
Bases:
simsoptpp.SurfaceRZFourier
,simsopt.geo.surface.Surface
SurfaceRZFourier
is a surface that is represented in cylindrical coordinates using the following Fourier series:\[r(\theta, \phi) = \sum_{m=0}^{m_{\text{pol}}} \sum_{n=-n_{\text{tor}}}^{n_\text{tor}} [ r_{c,m,n} \cos(m \theta - n_{\text{fp}} n \phi) + r_{s,m,n} \sin(m \theta - n_{\text{fp}} n \phi) ]\]and the same for \(z(\theta, \phi)\).
Here, \((r,\phi, z)\) are standard cylindrical coordinates, and theta is any poloidal angle.
Note that for \(m=0\) we skip the \(n<0\) term for the cos terms, and the \(n \leq 0\) for the sin terms.
In addition, in the
stellsym=True
case, we skip the sin terms for \(r\), and the cos terms for \(z\).For more information about the arguments
nphi
,ntheta
,range
,quadpoints_phi
, andquadpoints_theta
, see the general documentation on Surfaces.- Parameters
nfp – The number of field periods.
stellsym – Whether the surface is stellarator-symmetric, i.e. symmetry under rotation by \(\pi\) about the x-axis.
mpol – Maximum poloidal mode number included.
ntor – Maximum toroidal mode number included, divided by
nfp
.nphi – Number of grid points \(\phi_j\) in the toroidal angle \(\phi\).
ntheta – Number of grid points \(\theta_j\) in the toroidal angle \(\theta\).
range – Toroidal extent of the \(\phi\) grid. Set to
"full torus"
(or equivalentlySurfaceRZFourier.RANGE_FULL_TORUS
) to generate points up to 1 (with no point at 1). Set to"field period"
(or equivalentlySurfaceRZFourier.RANGE_FIELD_PERIOD
) to generate points up to \(1/n_{fp}\) (with no point at \(1/n_{fp}\)). Set to"half period"
(or equivalentlySurfaceRZFourier.RANGE_HALF_PERIOD
) to generate points up to \(1/(2 n_{fp})\) (with no point at \(1/(2 n_{fp})\)). Ifquadpoints_phi
is specified,range
is irrelevant.quadpoints_phi – Set this to a list or 1D array to set the \(\phi_j\) grid points directly.
quadpoints_theta – Set this to a list or 1D array to set the \(\theta_j\) grid points directly.
- _make_names()
Form a list of names of the rc, zs, rs, or zc array elements.
- _validate_mn(m, n)
Check whether m and n are in the allowed range.
- change_resolution(mpol, ntor)
Change the values of mpol and ntor. Any new Fourier amplitudes will have a magnitude of zero. Any previous nonzero Fourier amplitudes that are not within the new range will be discarded.
- darea()
Short hand for Surface.darea_by_dcoeff()
- dvolume()
Short hand for Surface.dvolume_by_dcoeff()
- fixed_range(mmin, mmax, nmin, nmax, fixed=True)
Set the ‘fixed’ property for a range of m and n values.
All modes with m in the interval [mmin, mmax] and n in the interval [nmin, nmax] will have their fixed property set to the value of the fixed parameter. Note that mmax and nmax are included (unlike the upper bound in python’s range(min, max).)
- classmethod from_focus(filename, **kwargs)
Read in a surface from a FOCUS-format file.
- classmethod from_vmec_input(filename: str, **kwargs)
Read in a surface from a VMEC input file. The
INDATA
namelist of this file will be read using f90nml. Note that this function does not require the VMEC python module.- Parameters
filename – Name of the
input.*
file to read.kwargs – Any other arguments to pass to the
SurfaceRZFourier
constructor. You can specifyquadpoints_theta
andquadpoints_phi
here.
- classmethod from_wout(filename: str, s: float = 1.0, interp_kind: str = 'linear', **kwargs)
Read in a surface from a VMEC wout output file. Note that this function does not require the VMEC python module.
- Parameters
filename – Name of the
wout_*.nc
file to read.s – Value of normalized toroidal flux to use for the surface. The default value of 1.0 corresponds to the VMEC plasma boundary. Must lie in the interval [0, 1].
interp_kind – Interpolation method in s. The available options correspond to the
kind
argument of scipy.interpolate.interp1d.kwargs – Any other arguments to pass to the
SurfaceRZFourier
constructor. You can specifyquadpoints_theta
andquadpoints_phi
here.
- get_dofs()
Return the dofs associated to this surface.
- get_rc(m, n)
Return a particular rc Parameter.
- get_rs(m, n)
Return a particular rs Parameter.
- get_zc(m, n)
Return a particular zc Parameter.
- get_zs(m, n)
Return a particular zs Parameter.
- recompute_bell(parent=None)
Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.
This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.
Need to be implemented by classes that provide a dof_setter for external handling of DOFs.
- return_fn_map: Dict[str, Callable] = {'area': <instancemethod area>, 'aspect-ratio': <function Surface.aspect_ratio>, 'volume': <instancemethod volume>}
- set_dofs(self: simsoptpp.SurfaceRZFourier, arg0: List[float]) None
- set_rc(m, n, val)
Set a particular rc Parameter.
- set_rs(m, n, val)
Set a particular rs Parameter.
- set_zc(m, n, val)
Set a particular zc Parameter.
- set_zs(m, n, val)
Set a particular zs Parameter.
- to_RZFourier()
No conversion necessary.
- write_nml(filename: str = 'boundary')
Writes a fortran namelist file containing the RBC/RBS/ZBS/ZBS coefficients, in the form used in VMEC and SPEC input files.
- Parameters
filename – Name of the file to write.
simsopt.geo.surfacexyzfourier module
- class simsopt.geo.surfacexyzfourier.SurfaceXYZFourier(nfp=1, stellsym=True, mpol=1, ntor=0, nphi=None, ntheta=None, range='full torus', quadpoints_phi=None, quadpoints_theta=None)
Bases:
simsoptpp.SurfaceXYZFourier
,simsopt.geo.surface.Surface
SurfaceXYZFourier is a surface that is represented in Cartesian coordinates using the following Fourier series:
\[\begin{split}\hat x(\phi,\theta) &= \sum_{m=0}^{m_\text{pol}} \sum_{n=-n_{\text{tor}}}^{n_{tor}} [ x_{c,m,n} \cos(m \theta - n_\text{ fp } n \phi) + x_{s,m,n} \sin(m \theta - n_\text{ fp } n \phi)]\\ \hat y(\phi,\theta) &= \sum_{m=0}^{m_\text{pol}} \sum_{n=-n_\text{tor}}^{n_\text{tor}} [ y_{c,m,n} \cos(m \theta - n_\text{fp} n \phi) + y_{s,m,n} \sin(m \theta - n_\text{fp} n \phi)]\\ z(\phi,\theta) &= \sum_{m=0}^{m_\text{pol}} \sum_{n=-n_\text{tor}}^{n_\text{tor}} [ z_{c,m,n} \cos(m \theta - n_\text{fp}n \phi) + z_{s,m,n} \sin(m \theta - n_\text{fp}n \phi)]\end{split}\]where
\[\begin{split}x &= \hat x \cos(\phi) - \hat y \sin(\phi)\\ y &= \hat x \sin(\phi) + \hat y \cos(\phi)\end{split}\]Note that for \(m=0\) we skip the \(n<0\) term for the cos terms, and the \(n \leq 0\) for the sin terms.
When enforcing stellarator symmetry, we set the
\[x_{s,*,*}, ~y_{c,*,*}, \text{and} ~z_{c,*,*}\]terms to zero.
For more information about the arguments
nphi
,ntheta
,range
,quadpoints_phi
, andquadpoints_theta
, see the general documentation on Surfaces.- Parameters
nfp – The number of field periods.
stellsym – Whether the surface is stellarator-symmetric, i.e. symmetry under rotation by \(\pi\) about the x-axis.
mpol – Maximum poloidal mode number included.
ntor – Maximum toroidal mode number included, divided by
nfp
.nphi – Number of grid points \(\phi_j\) in the toroidal angle \(\phi\).
ntheta – Number of grid points \(\theta_j\) in the toroidal angle \(\theta\).
range – Toroidal extent of the \(\phi\) grid. Set to
"full torus"
(or equivalentlySurfaceXYZFourier.RANGE_FULL_TORUS
) to generate points up to 1 (with no point at 1). Set to"field period"
(or equivalentlySurfaceXYZFourier.RANGE_FIELD_PERIOD
) to generate points up to \(1/n_{fp}\) (with no point at \(1/n_{fp}\)). Set to"half period"
(or equivalentlySurfaceXYZFourier.RANGE_HALF_PERIOD
) to generate points up to \(1/(2 n_{fp})\) (with no point at \(1/(2 n_{fp})\)). Ifquadpoints_phi
is specified,range
is irrelevant.quadpoints_phi – Set this to a list or 1D array to set the \(\phi_j\) grid points directly.
quadpoints_theta – Set this to a list or 1D array to set the \(\theta_j\) grid points directly.
- get_dofs()
Return the dofs associated to this surface.
- recompute_bell(parent=None)
Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.
This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.
Need to be implemented by classes that provide a dof_setter for external handling of DOFs.
- return_fn_map: Dict[str, Callable] = {'area': <instancemethod area>, 'aspect-ratio': <function Surface.aspect_ratio>, 'volume': <instancemethod volume>}
- set_dofs(dofs)
Set the dofs associated to this surface.
- to_RZFourier()
Return a SurfaceRZFourier instance corresponding to the shape of this surface.
simsopt.geo.surfacexyztensorfourier module
- class simsopt.geo.surfacexyztensorfourier.SurfaceXYZTensorFourier(nfp=1, stellsym=True, mpol=1, ntor=1, clamped_dims=[False, False, False], nphi=None, ntheta=None, range='full torus', quadpoints_phi=None, quadpoints_theta=None)
Bases:
simsoptpp.SurfaceXYZTensorFourier
,simsopt.geo.surface.Surface
SurfaceXYZTensorFourier is a surface that is represented in cartesian coordinates using the following Fourier series:
\[\begin{split}\hat x(\theta, \phi) &= \sum_{i=0}^{2m_\text{pol}} \sum_{j=0}^{2n_\text{tor}} x_{ij} w_i(\theta)v_j(\phi)\\ \hat y(\theta, \phi) &= \sum_{i=0}^{2m_\text{pol}} \sum_{j=0}^{2n_\text{tor}} y_{ij} w_i(\theta)v_j(\phi)\\ x(\phi, \theta) &= \hat x(\phi, \theta) \cos(\phi) - \hat y(\phi, \theta) \sin(\phi)\\ y(\phi, \theta) &= \hat x(\phi, \theta) \sin(\phi) + \hat y(\phi, \theta) \cos(\phi)\\ z(\theta, \phi) &= \sum_{i=0}^{2m_\text{pol}} \sum_{j=0}^{2n_\text{tor}} z_{ij} w_i(\theta)v_j(\phi)\end{split}\]where the basis functions \({v_j}\) are given by
\[\{1, \cos(1\,\mathrm{nfp}\,\phi), \ldots, \cos(n_\text{tor}\,\mathrm{nfp}\,\phi), \sin(1\,\mathrm{nfp}\,\phi), \ldots, \sin(n_\text{tor}\,\mathrm{nfp}\,\phi)\}\]and \({w_i}\) are given by
\[\{1, \cos(1\theta), \ldots, \cos(m_\text{pol}\theta), \sin(1\theta), \ldots, \sin(m_\text{pol}\theta)\}\]When stellsym=True the sums above change as follows:
\[\begin{split}\hat x(\theta, \phi) &= \sum_{i=0}^{m_\text{pol}} \sum_{j=0}^{n_\text{tor}} x_{ij} w_i(\theta)v_j(\phi) + \sum_{i=m_\text{pol}+1}^{2m_\text{pol}} \sum_{j=n_\text{tor}+1}^{2n_\text{tor}} x_{ij} w_i(\theta)v_j(\phi)\\ \hat y(\theta, \phi) &= \sum_{i=0}^{m_\text{pol}} \sum_{j=n_\text{tor}+1}^{2n_\text{tor}} y_{ij} w_i(\theta)v_j(\phi) + \sum_{i=m_\text{pol}+1}^{2m_\text{pol}} \sum_{j=0}^{n_\text{tor}} y_{ij} w_i(\theta)v_j(\phi)\\\\ z(\theta, \phi) &= \sum_{i=0}^{m_\text{pol}} \sum_{j=n_\text{tor}+1}^{2n_\text{tor}} z_{ij} w_i(\theta)v_j(\phi) + \sum_{i=m_\text{pol}+1}^{2m_\text{pol}} \sum_{j=0}^{n_\text{tor}} z_{ij} w_i(\theta)v_j(\phi)\end{split}\]For more information about the arguments
nphi
,ntheta
,range
,quadpoints_phi
, andquadpoints_theta
, see the general documentation on Surfaces.- Parameters
nfp – The number of field periods.
stellsym – Whether the surface is stellarator-symmetric, i.e. symmetry under rotation by \(\pi\) about the x-axis.
mpol – Maximum poloidal mode number included.
ntor – Maximum toroidal mode number included, divided by
nfp
.clamped_dims –
???
nphi – Number of grid points \(\phi_j\) in the toroidal angle \(\phi\).
ntheta – Number of grid points \(\theta_j\) in the toroidal angle \(\theta\).
range – Toroidal extent of the \(\phi\) grid. Set to
"full torus"
(or equivalentlySurfaceXYZTensorFourier.RANGE_FULL_TORUS
) to generate points up to 1 (with no point at 1). Set to"field period"
(or equivalentlySurfaceXYZTensorFourier.RANGE_FIELD_PERIOD
) to generate points up to \(1/n_{fp}\) (with no point at \(1/n_{fp}\)). Set to"half period"
(or equivalentlySurfaceXYZTensorFourier.RANGE_HALF_PERIOD
) to generate points up to \(1/(2 n_{fp})\) (with no point at \(1/(2 n_{fp})\)). Ifquadpoints_phi
is specified,range
is irrelevant.quadpoints_phi – Set this to a list or 1D array to set the \(\phi_j\) grid points directly.
quadpoints_theta – Set this to a list or 1D array to set the \(\theta_j\) grid points directly.
- get_dofs()
Return the dofs associated to this surface.
- get_stellsym_mask()
In the case of stellarator symmetry, some of the information is redundant, since the coordinates at (-phi, -theta) are the same (up to sign changes) to those at (phi, theta). The point of this function is to identify those angles phi and theta that we can ignore. This is difficult to do in general, so we focus on the following three common cases below:
phis = np.linspace(0, 1/self.nfp, 2*ntor+1, endpoint=False) thetas = np.linspace(0, 1, 2*mpol+1, endpoint=False)
or
phis = np.linspace(0, 1/self.nfp, 2*ntor+1, endpoint=False) thetas = np.linspace(0, 0.5, 2*mpol, endpoint=False)
or
phis = np.linspace(0, 1/(2*self.nfp), ntor+1, endpoint=False) thetas = np.linspace(0, 1, 2*mpol+1, endpoint=False)
This function could be extended to be aware of rotational symmetry as well. So far we assume that that redundancy was removed already (hence the phis only go to 1/nfp or 1/(2*nfp)).
- recompute_bell(parent=None)
Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.
This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.
Need to be implemented by classes that provide a dof_setter for external handling of DOFs.
- set_dofs(dofs)
Set the dofs associated to this surface.
- to_RZFourier()
Return a SurfaceRZFourier instance corresponding to the shape of this surface.