Adjoint approach to calculating shape gradients for 3D magnetic confinement equilibria

TL;DR

This paper introduces an efficient adjoint method for calculating shape gradients in 3D magnetic confinement equilibria, significantly reducing computational effort for optimization and sensitivity analysis.

Contribution

It presents a novel adjoint approach leveraging self-adjoint properties of MHD equations to compute shape gradients more efficiently in 3D toroidal equilibria.

Findings

Reduces the number of equilibrium computations by factors of O(N)

Applicable to shape changes of flux surfaces and coils

Validated with multiple figures of merit

Abstract

The shape gradient quantifies the change in some figure of merit resulting from differential perturbations to a shape. Shape gradients can be applied to gradient-based optimization, sensitivity analysis, and tolerance calculation. An efficient method for computing the shape gradient for toroidal 3D MHD equilibria is presented. The method is based on the self-adjoint property of the equations for driven perturbations of MHD equilibria and is similar to the Onsager symmetry of transport coefficients. Two versions of the shape gradient are considered. One describes the change in a figure of merit due to an arbitrary displacement of the outer flux surface; the other describes the change in the figure of merit due to the displacement of a coil. The method is implemented for several example figures of merit and compared with direct calculation of the shape gradient. In these examples the adjoint method reduces the number of equilibrium computations by factors of $\mathcal{O}(N)$, where $N$ is the number of parameters used to describe the outer flux surface or coil shapes.

Figures

Equations

$$\delta f(S_P;\delta \text{r})={\int }_{S_P}{d}^{2}xS\delta \text{r}\cdot \text{n}.$$

$$\delta f(C;\delta {\text{r}}_C)=\sum _k{\int }_{C_k}dl{\text{S}}_k\cdot \delta {\text{r}}_{C_k}.$$

$$\text{B}=\nabla \psi \times \nabla \theta -\nabla \Phi \times \nabla \zeta =\nabla \psi \times \nabla \alpha .$$

$$\alpha =\theta -\iota (\psi )\zeta ,$$

$$0=-\nabla p(\psi )+\frac{\text{J}\times \text{B}}{c},$$

$$\nabla \times \text{B}=\frac{4\pi }{c}\text{J},$$

$$I_T(\psi )={\int }_{S_T(\psi )}{d}^{2}x\text{n}\cdot \text{J}={\int }_{S_T(\psi )}d\psi d\theta \sqrt{g}\nabla \zeta \cdot \text{J},$$

$$\delta \text{B}=\nabla \delta \psi \times \nabla \theta +\nabla \psi \times \nabla \delta \theta -\nabla \Phi \times \nabla \delta \zeta -\nabla \left(\iota \delta \psi +\delta \Phi \right)\times \nabla \zeta .$$

$$\delta \text{B}=\nabla \times \left(\delta \psi \nabla \theta -\iota \delta \psi \nabla \zeta -\delta \theta \nabla \psi +\delta \zeta \nabla \Phi \right)-\nabla \delta \Phi (\psi )\times \nabla \zeta .$$

$$\delta \text{A}=\mathbf{\xi }\times \text{B}-\delta \Phi \nabla \zeta .$$

$$\mathbf{\xi }\times \text{B}=\delta \psi \left(\nabla \theta -\iota \nabla \zeta \right)-\delta \theta \nabla \psi +\delta \zeta \nabla \Phi ,$$

$$\delta I_T(\psi )={\int }_{\partial S_T(\psi )}d\theta \sqrt{g}\mathbf{\xi }\cdot \nabla \psi \text{J}\cdot \nabla \zeta +{\int }_{S_T(\psi )}d\psi d\theta \sqrt{g}\delta \text{J}\cdot \nabla \zeta ,$$

$$U_T=\frac{1}{c}{\int }_{V_T}{d}^{3}x\left(\delta {\text{J}}_1\cdot \delta {\text{A}}_2-\delta {\text{J}}_2\cdot \delta {\text{A}}_1\right)=0,$$

$$U_T=U_P+U_B+U_C=0.$$

$$U_P=\frac{1}{c}{\int }_{V_P}{d}^{3}x\left(\delta {\text{J}}_1\cdot \left({\mathbf{\xi }}_2\times \text{B}-\delta {\Phi }_2\nabla \zeta \right)-\delta {\text{J}}_2\cdot \left({\mathbf{\xi }}_1\times \text{B}-\delta {\Phi }_1\nabla \zeta \right)\right).$$

$$0=\delta \text{F}+\nabla (\mathbf{\xi }\cdot \nabla p)+\frac{\delta \text{J}\times \text{B}+\text{J}\times \delta \text{B}}{c},$$

$$U_B=\frac{1}{c}{\int }_{S_p}\frac{d^{2}x}{|\nabla \psi |}\left({\mathbf{\xi }}_1\cdot \nabla \psi \text{J}\cdot \delta {\text{A}}_2-{\mathbf{\xi }}_2\cdot \nabla \psi \text{J}\cdot \delta {\text{A}}_1\right),$$

$$U_B=\frac{1}{c}{\int }_{S_p}d\theta d\zeta \sqrt{g}\text{J}\cdot \nabla \zeta \left(-{\mathbf{\xi }}_1\cdot \nabla \psi \delta {\Phi }_2+{\mathbf{\xi }}_2\cdot \nabla \psi \delta {\Phi }_1\right).$$

$$U_C=\frac{1}{c}{\int }_{V_V}{d}^{3}x\left(\delta {\text{J}}_{C_1}\cdot \delta {\text{A}}_{V_2}-\delta {\text{J}}_{C_2}\cdot \delta {\text{A}}_{V_1}\right),$$

$$\begin{gathered} {\int }_{V_P}{d}^{3}x\left(-{\mathbf{\xi }}_1\cdot \delta {\text{F}}_2+{\mathbf{\xi }}_2\cdot \delta {\text{F}}_1\right)+\frac{2\pi }{c}{\int }_{V_P}d\psi \left(\delta {\Phi }_1\frac{d\delta I_{T,2}}{d\psi }-\delta {\Phi }_2\frac{d\delta I_{T,1}}{d\psi }\right) \\ +\frac{1}{c}{\int }_{V_V}{d}^{3}x\left(\delta {\text{J}}_{C_1}\cdot \delta {\text{A}}_{V_2}-\delta {\text{J}}_{C_2}\cdot \delta {\text{A}}_{V_1}\right)=0. \end{gathered}$$

$$U_P=\frac{1}{4\pi }{\int }_{S_P}{d}^{2}x\text{n}\cdot \left(\delta {\text{B}}_1\times \delta {\text{A}}_2-\delta {\text{B}}_2\times \delta {\text{A}}_1\right).$$

$$\begin{gathered} {\int }_{V_P}{d}^{3}x\left(-{\mathbf{\xi }}_1\cdot \delta {\text{F}}_2+{\mathbf{\xi }}_2\cdot \delta {\text{F}}_1\right)-\frac{2\pi }{c}{\int }_{V_P}d\psi \left(\delta I_{T,2}\frac{d\delta {\Phi }_1}{d\psi }-\delta I_{T,1}\frac{d\delta {\Phi }_2}{d\psi }\right) \\ -\frac{1}{4\pi }{\int }_{S_P}{d}^{2}x\text{n}\cdot \left({\mathbf{\xi }}_2\delta {\text{B}}_1\cdot \text{B}-{\mathbf{\xi }}_1\delta {\text{B}}_2\cdot \text{B}\right)=0. \end{gathered}$$

$$\frac{2\pi }{c}{\int }_{V_P}d\psi \left(\delta I_{T,2}\frac{d\delta {\Phi }_1}{d\psi }-\delta I_{T,1}\frac{d\delta {\Phi }_2}{d\psi }\right)=\frac{1}{4\pi }{\int }_{S_p}{d}^{2}x\text{n}\cdot {\mathbf{\xi }}_1\delta {\text{B}}_2\cdot \text{B}.$$

$$\frac{2\pi }{c}{\int }_{V_P}d\psi \left(\delta I_{T,2}\frac{d\delta {\Phi }_1}{d\psi }\right)=\frac{1}{4\pi }{\int }_{S_P}{d}^{2}x\mathbf{n}\cdot {\mathbf{\xi }}_1\delta {\text{B}}_2\cdot \text{B}.$$

$$\frac{2\pi }{c}{\int }_{V_P}d\psi \left(\delta I_{T,2}\frac{d\delta {\Phi }_1}{d\psi }\right)=\frac{1}{c}{\int }_{V_V}{d}^{3}x\delta {\text{J}}_{C_1}\cdot \delta {\text{A}}_{V_2},$$

$$\frac{1}{c}{\int }_{V_V}{d}^{3}x\delta {\text{J}}_{C_1}\cdot \delta {\text{A}}_{V_2}=\frac{1}{c}\sum _k\left(\delta {\Phi }_{C_{2,k}}\delta I_{C_{1,k}}+I_{C_k}{\int }_{C_k}dl\delta {\text{r}}_{C_k}(\text{x})\cdot \text{t}\times \delta {\text{B}}_2\right).$$

$${\int }_{V_P}{d}^{3}x{\mathbf{\xi }}_1\cdot \delta {\text{F}}_2=\frac{1}{c}{\int }_{V_V}{d}^{3}x\delta {\text{J}}_{C_1}\cdot \delta {\text{A}}_{V_2}.$$

$$f_{\beta }=\frac{f_P}{f_B},$$

$$f_P={\int }_{V_p}{d}^{3}xp(\psi )$$

$$f_B={\int }_{V_p}{d}^{3}x\frac{B^2}{8\pi }.$$

Full text

1 Introduction

The design of stellarator magnetohydrodynamic (MHD) equilibria requires optimizing within a high dimensional space due to the fully 3-dimensional nature of the magnetic field configuration and the sensitive dependence of charged particle trajectories in such field configurations (Boozer, 2015). While there are many possible choices for the space in which to optimize, a common choice is the space of the shape of the outer boundary of the plasma (Hirshman et al., 1999; Drevlak et al., 2018). The confining electromagnetic coils must then be designed to reproduce the desired plasma boundary. This approach was used to design Wendelstein 7-X (Grieger et al., 1992) and the Helically Symmetric Experiment (HSX) (Anderson et al., 1995). An alternative approach is to optimize the shape of electromagnetic coils directly to minimize an objective function which includes functions of the equilibria (Hanson & Cary, 1984). This can be performed with the merged STELLOPT/COILOPT code, which was used in the late stages of the NCSX design (Strickler et al., 2004). To navigate such spaces, it is often useful to employ gradient-based optimization techniques. For this reason, it is desirable to compute derivatives with respect to shape.

As the target optimized configuration can never be realized exactly, an analysis of the sensitivity to perturbations, such as errors in coil fabrication or assembly, is central to the success of a stellarator. Tight tolerances have proven to be a significant driver of the cost of stellarator experiments (Strykowsky et al., 2009; Klinger et al., 2013); thus an improvement to the algorithms used to conduct sensitivity studies can have great impact on the field. To quantify the coil tolerances for flux surface quality of LHD (Yamazaki et al., 1993) and NCSX (Brooks & Reiersen, 2003; Williamson et al., 2005), perturbations of several distributions were manually applied to each coil. Sensitivity analysis can also be performed with analytic derivatives. Numerical derivatives with respect to tilt angle and coil translation of the CNT coils have been used to compute the sensitivity of the rotational transform on axis (Hammond et al., 2016). Analytic derivatives have recently been applied to study coil sensitivities of the CNT stellarator by considering the eigenvectors of the Hessian matrix (Zhu et al., 2018). Thus, in addition to gradient-based optimization, derivatives with respect to shape can be applied to sensitivity analysis.

The gradients of figures of merit with respect to shape have often been represented as derivatives with respect to whichever quantities parameterize the shape. Examples of such quantities are the amplitudes ($R_{mn}^{c}$, $Z_{mn}^{s}$) in the double Fourier series for the cylindrical coordinates of a toroidal surface. Another way to represent derivatives with respect to shape is the shape gradient (Landreman & Paul, 2018), which provides a local and coordinate-independent form. Consider any scalar figure of merit, $f$, that depends on a 3D MHD equilibrium solution. We can consider $f$ to be a functional of the shape of the outer boundary of the plasma, $S_{P}$. In this case, a differential change to the boundary, $\delta\textbf{r}$, causes a corresponding change to the figure of merit, $\delta f$,

[TABLE]

Here n is the outward unit normal and $S$ is the shape gradient. The shape gradient quantifies the local linear sensitivity of a figure of merit to differential perturbations of the shape. As tangential displacements to $S_{P}$ do not cause any changes to $f$, $\delta f$ only depends on the normal component of $\delta\textbf{r}$. The Hadamard-Zolésio structure theorem (Delfour & Zolésio, 2011) states that under certain assumptions of smoothness, the shape derivative of a functional can be written in the form of (1). This can be thought of an instance of the Riesz representation theorem, which states that any linear functional, such as $\delta f(\delta\textbf{r})$, can be written as an inner product over the appropriate space (Rudin, 2006). The motivation for the form of (1) is discussed in more detail in section 2 of Landreman & Paul (2018).

If we consider $f$ to be a functional of the shape of the electromagnetic coils, a differential change to the coils, $\delta\textbf{r}_{C}$, will cause a corresponding change to the figure of merit, $\delta f$,

[TABLE]

Here the sum is taken over the coils, $C_{k}$ is a curve describing the filamentary shape of each coil, $C=\{C_{k}\}$, and $\delta\textbf{r}_{C}=\{\delta\textbf{r}_{C_{k}}\}$. The shape gradient for coil $k$ is $\textbf{S}_{k}$. As tangential displacements to the coil do not change $f$, $\textbf{S}_{k}$ must be perpendicular to the tangent vector along $C_{k}$.

As the shape gradient provides local sensitivity information, it can be used to quantify engineering tolerances with respect to the displacement of coils or magnetic perturbations. The shape gradient representation can be computed from parameter derivatives by solving a small linear system (Landreman & Paul, 2018).

However, computing parameter derivatives can often be computationally expensive, as numerical derivatives require evaluating the objective function at least $N+1$ times for $N$ parameters if one-sided finite differences are used, or $2N$ times for centered differences. As computing the objective function often involves solving a linear or nonlinear system, such as the MHD equilibrium equations, this implies solving the system of equations $\geq N+1$ times. Numerical derivatives also introduce additional noise, and the finite difference step size must be chosen carefully. To avoid these difficulties, it is advantageous to compute shape gradients using adjoint methods (Pironneau, 1974; Glowinski & Pironneau, 1975; Dekeyser et al., 2012, 2014a, 2014b). Adjoint methods allow the analytic derivative with respect to all $N$ parameters to be computed with only two solutions to the system of equations. Adjoint methods are thus much more efficient for computing derivatives with respect to many parameters, and they do not introduce the noise of numerical derivatives. Adjoint methods were recently used in the context of stellarator design by Paul et al. (2018) for shape optimization of coil winding surfaces.

Rather than use parameter derivatives, in this work we will use an adjoint method to compute the shape gradient directly. This is sometimes termed adjoint shape sensitivity or adjoint shape optimization, which has its origins in aerodynamic engineering and computational fluid dynamics (Pironneau, 1974; Glowinski & Pironneau, 1975). As with adjoint methods for parameter derivatives, this technique only requires the solution of two linear or non-linear systems of equations. This technique has been applied to magnetic confinement fusion for the design of tokamak divertor shapes by solving forward and adjoint fluid edge equations (Dekeyser et al., 2012, 2014a, 2014b). As stellarators require many parameters to fully describe their shape, adjoint shape sensitivity could significantly decrease the cost of computing the shape gradient. If one is optimizing in the space of parameters describing the boundary of the plasma or the shape of coils, the shape gradient representation obtained from the adjoint method can be converted to parameter derivatives upon multiplication with a small matrix (Landreman & Paul, 2018).

In section 2, the fundamental adjoint relations for perturbations to MHD equilibria are derived and discussed. These relations take a form that is similar to that of transport coefficients that are related by Onsager symmetry (Onsager, 1931). Specifically, perturbations to the equilibrium are characterized as a set of generalized responses to a complementary set of generalized forces. The responses and forces can be thought of as being related by a matrix operator, which is symmetric. The resulting relations among forces and responses can be used to compute the shape gradient of functions of the equilibria with respect to displacements of the plasma boundary, as in 1, or the coil shapes, as in 2, Several applications to stellarator figures of merit will be demonstrated in section 3. While the primary application considered in this work will be stellarator optimization, the relations we obtain are equally applicable for 2D equilibria.

2 Driven linear perturbations of 3D MHD equilibria

The goal is to find a relation between small perturbations to a 3D magnetic confinement equilibrium configuration and the resulting change in a figure of merit of interest. As mentioned, the perturbations may be prescribed in one of two ways: either as a change in the shape of the outermost flux surface (fixed boundary) or as a change in the position, shape, or current strength in the coils confining the plasma (free boundary). (Even though the boundary shape changes in the former case, we refer to it as “fixed boundary” since the equilibrium code is run in fixed boundary mode, and since the associated adjoint problem will turn out to have no boundary perturbation.) These perturbations will be referred to as the “true” perturbations and their direct calculation needs to be repeated many times for each possible change in shape of the outer flux surface or each change in the coil configuration to determine fully the sensitivity of the equilibrium.

The approach we use is to instead calculate a different change in the equilibrium, which we refer to as the “adjoint” perturbation. The adjoint perturbation will correspond to the change in the equilibrium when an additional bulk force acts on the plasma, or the toroidal current profile is changed. For the adjoint perturbation there is no change to the outer flux surface in the fixed boundary case or to the coil currents in the free boundary case. In this section we will show that aspects of the true and adjoint changes are related to each other in a manner similar to Onsager symmetry. Thus, it will be shown that by calculating the adjoint perturbation, with a judiciously chosen added force or change in the toroidal current profile, the solution to the true problem can be determined.

We consider equilibria in which the magnetic field in the plasma can be expressed in terms of scalar of functions $\psi(\textbf{x}),\Phi(\psi),\theta(\textbf{x})$, and $\zeta(\textbf{x})$,

[TABLE]

We will regard $\psi$ as labeling the flux surfaces and consider toroidal geometries for which

[TABLE]

label field lines in a flux surface, where $\theta$ is a poloidal angle, $\zeta$ is a toroidal angle, and $\iota(\psi)=d\Phi/d\psi$ is the rotational transform, with $\Phi$ being the poloidal flux function. (Any straight-field-line angles may be used.) With these definitions the magnetic flux passing toroidally through a poloidally closed curve of constant $\psi$ is $2\pi\psi$. The flux passing poloidally between the magnetic axis and the surface of constant $\psi$ is $2\pi\Phi(\psi)$. Thus, we assume that good flux surfaces exist and leave aside the issues of islands and chaotic field lines.

In addition to the representation of the magnetic field, we assume that MHD force balance is satisfied with a scalar pressure, $p(\psi)$,

[TABLE]

where the current density, J, satisfies Ampere’s law,

[TABLE]

and $c$ is the speed of light. As mentioned, we will consider two cases, a fixed boundary case in which the shape of the outer flux surface is prescribed, and a free boundary case for which outside the plasma, whose surface is defined by a particular value of toroidal flux, the force balance equation (5) does not apply, but rather, the magnetic field is determined by Ampere’s law (6) with a given current density $\textbf{J}_{c}$, representing current flowing in a set of coils.

From (5) it follows that current density stream-lines also lie in the $\psi=$ constant surfaces. The toroidal current passing through a surface, $S_{T}(\psi)$, whose perimeter is a closed poloidal loop at constant $\psi$ is given by

[TABLE]

where $\sqrt{g}^{-1}=\nabla\zeta\cdot\nabla\psi\times\nabla\theta$.

Equations 3, 4, 4, 5, 6 and 7 describe our base equilibrium configuration. We now consider small changes in the equilibrium that are assumed to yield a second equilibrium state of the same form as 3, but with new functions such that $\textbf{B}^{\prime}=\nabla\psi^{\prime}\times\nabla\theta^{\prime}-\nabla\Phi^{\prime}(\psi^{\prime})\times\nabla\zeta^{\prime}$. Each of the primed variables is assumed to differ from the corresponding unprimed variables by a small amount (e.g. $\psi^{\prime}=\psi+\delta\psi(\textbf{x})$). The perturbed magnetic field can then be expressed $\textbf{B}^{\prime}=\textbf{B}+\delta\textbf{B}$, where

[TABLE]

We write the perturbed poloidal flux as the sum of a term resulting from the perturbation of toroidal flux at fixed rotational transform, $\iota\delta\psi$, and a term representing the perturbed rotational transform, $\delta\Phi(\psi)$. Thus, we can regroup the terms in 8 as follows

[TABLE]

The group of terms in parentheses in 9 corresponds to perturbations of the magnetic field allowed by ideal MHD, which is constrained by the “frozen in law”, and which preserves the rotational transform, ($\delta\iota=0$). The last term in 9 allows for changes in the rotational transform, ($\delta\iota=d\delta\Phi/d\psi$). Note also that the expression in parentheses in 9 can be written as a sum of terms parallel to $\nabla\psi$ and $\nabla\alpha$, and hence it is perpendicular to B. The group of terms in parentheses in 9 can thus be expressed in terms of a vector potential that is perpendicular to the equilibrium magnetic field, while the last term in 9 can be represented in terms of a vector potential in the toroidal direction, which thus has a component parallel to the equilibrium field. We can therefore write $\delta\textbf{B}=\nabla\times\delta\textbf{A}$, where

[TABLE]

Here, the variable $\bm{\xi}$ can be taken to be perpendicular to the applied magnetic field, and can be thought of as a displacement of the equilibrium magnetic field line. Using 8 we write

[TABLE]

and from this we can see that perturbations of the toroidal and poloidal angles correspond to displacements in the flux surface, and the perturbation $\delta\psi$ gives a displacement out of the flux surface. A surface containing a prescribed toroidal flux $2\pi\psi_{0}$ in the unperturbed field is defined by $\psi(\textbf{x})=\psi_{0}$. From 11, the perturbation to $\psi$ at fixed position is $\delta\psi=-\bm{\xi}\cdot\nabla\psi$. When the field is perturbed the surface moves. As $\bm{\xi}$ is a vector which describes the motion of field lines, the perturbed toroidal flux label moving with the $\psi$ surface, as measured in the unperturbed coordinate system, is given by $\psi(\textbf{x})=\psi_{0}+\bm{\xi}\cdot\nabla\psi(\textbf{x})$.

The change in toroidal current flowing through the perturbed surface is

[TABLE]

where $S_{T}(\psi)$ is a surface at constant toroidal angle bounded by the $\psi$ surface and $\partial S_{T}(\psi)$ is the boundary of such surface, a closed poloidal loop. Here the first term accounts for the displacement of the flux surface and the second term accounts for the change in toroidal current density.

We now consider two distinct perturbations of the equilibrium of the type described by 10, 11 and 12, which we denote with subscripts 1 and 2. In general, variables with subscript 1 will be associated with the true perturbation, and those with subscripts 2 will be associated with the adjoint perturbation. We then form the quantity

[TABLE]

where the variables $\delta\textbf{J}$ and $\delta\textbf{A}$ are the changes in current density and vector potential associated with the two perturbations, and the integral is, for the time being, over all space. That $U_{T}=0$ follows from expressing the change in current densities in terms of the change in magnetic fields via Ampere’s law. This turns the integrand in 13 into a divergence, which in turn becomes a surface integral, which we take to be at infinity where fields vanish sufficiently fast.

We now express the volume integral in 13 as the sum of three terms,

[TABLE]

Here $U_{P}$ is the contribution from volume in the plasma, integrated just up to the plasma-vacuum boundary. For this term we represent the vector potentials using 10

[TABLE]

To evaluate 15 we use the perturbed force balance relation

[TABLE]

where $\delta\textbf{F}$ is an additional perturbed force to be prescribed. The term $\nabla(\bm{\xi}\cdot\nabla p)$ represents the perturbed force associated with the perturbed pressure in the case of a true solution for which the dependence of the pressure on flux is preserved. As the perturbation to the flux label at fixed position is $\delta\psi=-\bm{\xi}\cdot\nabla\psi$ and the pressure is assumed to be a fixed function $p(\psi)$, then the perturbation to the pressure at fixed position is $\delta p=-\bm{\xi}\cdot\nabla p(\psi)$. If the equilibrium calculation is performed with specified $p(\psi)$, this does not imply any additional restrictions.

The term $U_{B}$ comes from integrating over a thin layer at the plasma-vacuum boundary. At the boundary the difference between the perturbed and unperturbed current density has the character of a current sheet due to the displacement of the outermost flux surface. This effective current sheet causes a jump in the tangential components of the perturbation to the magnetic fields at the surface. This jump implies that care must be taken in evaluating the perturbed magnetic fields at the surface as they have different values on either side of the plasma-vacuum surface. However, the vector potential is continuous at the plasma-vacuum boundary. Thus, we write

[TABLE]

where the vector potentials are expressed as in 10. Using this expression for the vector potentials and expressing the surface integral as an integral over the toroidal and poloidal angles gives

[TABLE]

Here we note the terms in the vector potential coming from the MHD displacement cancel.

Last, the quantity $U_{C}$ represents the contribution from the integral over the volume outside the plasma where only the coil currents need to be included

[TABLE]

where $\delta\textbf{A}_{V_{1,2}}$ is the change in the vacuum vector potential, and $\delta\textbf{J}_{C_{1,2}}$ is the change in the coil current density.

Combining $U_{P}$, $U_{B}$, and $U_{C}$ gives the following relation appropriate to the free boundary case $U_{T}=U_{P}+U_{B}+U_{C}=0$, or

[TABLE]

The many steps leading to 20 are outlined in Appendix A.

A similar relation can be obtained in the fixed boundary case. Here the integral over the plasma volume, (15), can be written as a surface integral by applying the divergence theorem,

[TABLE]

Again, following steps outlined in Appendix A, this may be rewritten in the following form,

[TABLE]

The fixed boundary adjoint relation can also be obtained by applying the self-adjointness of the MHD force operator (Appendix B). If the second term in (22) is integrated by parts in $\psi$, we see that the fixed and free boundary adjoint relations share the terms involving the products of displacements with bulk forces and perturbed fluxes with perturbed toroidal currents. The integral over the vacuum region in (20) is replaced by an integral over the plasma boundary and a boundary term from the integration by parts in $\psi$ in (22).

We now have two integral relations between perturbations 1 and 2, viz 20 and 22. They have a common form in that they each are the sum of three integrals: the first involving forces and displacements, the second involving the toroidal current and poloidal flux profiles, and the third involving the manner in which the plasma boundary is prescribed. In 20, the free boundary case, the changes in coil current densities are specified. In 22, the fixed boundary case, the displacement of the outer flux surface is prescribed. Equations 20 and 22 can also be viewed as the difference in sums of generalized forces and responses. For example, in 20 we can consider the quantities $\delta\textbf{F}$, $\delta\Phi$, $\delta\textbf{J}_{c}$ as forces and $\bm{\xi}$, $d\delta I_{T}/d\psi$, $\delta\textbf{A}_{V}$ as responses. The fact that the sum of the products of true forces and adjoint responses less the products of adjoint forces and true responses vanishes is similar to the relation between forces and fluxes related by Onsager symmetry. In the case of Onsager symmetry this relation follows from the self-adjoint property of the collision operator. In the case of MHD equilibria it is known that the force operator is self-adjoint. Here we see that the self-adjoint property is extended to cases in which the MHD frozen-in constraint is broken.

These relations 20 and 22 can be used to generate the shape gradient if we manipulate the terms to be zero, or some known quantity. For example, if we impose that both the true and adjoint solutions maintain pressure as a function of flux, or equivalently there are no added forces, then the first terms in both 20 and 22 vanish. If we impose that the adjoint solution in 22 involves no change in shape of the outer flux surface, the boundary term involving $\textbf{n}\cdot\bm{\xi}_{2}$ vanishes. This leaves

[TABLE]

Thus, if we solve the adjoint problem with a prescribed change in toroidal current profile, and ask what is the change in rotational transform in the true problem when the current profile is unchanged we find

[TABLE]

The right-hand side has the form of 1, with $\delta\textbf{B}_{2}\cdot\textbf{B}/4\pi$ playing the role of the shape gradient. Thus, by picking an adjoint current profile change that is localized to a particular flux surface, we determine the sensitivity of the rotational transform at that flux surface to changes in the boundary shape. Note that even if the rotational transform is allowed to vary, the normal component of $\bm{\xi}_{1}$ describes the displacement of the boundary, as discussed in Appendix C. Thus the form of the shape gradient in (1) can be used with $\delta\textbf{r}\cdot\textbf{n}$ replaced with $\bm{\xi}_{1}\cdot\textbf{n}$.

If we impose that the adjoint solution in (20) involves no change in coil currents, upon integrating the middle term in 20 by parts a similar relation can be obtained for the free boundary case,

[TABLE]

where it has been assumed that $\delta I_{T,2}=0$ at the boundary.

When the coil currents are confined to wires, the right hand side can be expressed as changes in currents and fluxes and integrals along wires,

[TABLE]

Here $\delta\Phi_{C_{2,k}}$, $\delta I_{C_{1,k}}$ are the change in adjoint flux through and change in true current in coil $k$, and $I_{C_{k}}$ is the current through the unperturbed coil. The unit tangent vector along $C_{k}$ is t. The second term describes the effect of moving the coil. It is proportional to the current in the coil and the line integral of the displacement of the coil dotted with the cross product of the tangent to the coil and the perturbed adjoint magnetic field.

The sensitivity of other figures of merit may be evaluated by considering cases for which both the true and adjoint solutions preserve the rotational transform or the toroidal current profile. In this case a perturbed force density is included in the adjoint calculation

[TABLE]

In this case one must find an expression for the perturbed adjoint force, $\delta\textbf{F}_{2}$, such that the left hand side is the change in a figure of merit of interest. Some examples will be discussed in the next section.

3 Applications

In this section we will consider figures of merit which depend on the shape of the outer boundary of the plasma (section 3.1 and 3.2) and on the shape of the electromagnetic coils (section 3.3). The shape gradients of these figures of merit will be computed using both a direct method and an adjoint method, to demonstrate that the adjoint method produces identical results to the direct method but at much lower computational expense.

3.1 Surface shape gradient for $\beta$

Consider a figure of merit, the volume-averaged $\beta$,

[TABLE]

where

[TABLE]

and

[TABLE]

(This definition of volume-averaged $\beta$ is the one employed in the VMEC code (Hirshman & Whitson, 1983).) While $f_{\beta}$ is a figure of merit not often considered in stellarator shape optimization, we include this calculation to demonstrate the adjoint approach, as its shape gradient can be computed without modifications to an equilibrium code.

The differential change in $f_{P}$ associated with displacement $\bm{\xi}_{1}$ is

[TABLE]

The first term accounts for the change in $p$ at fixed position due to the motion of the flux surfaces, and the second term accounts for the motion of the boundary. The differential change in $f_{B}$ associated with $\bm{\xi}_{1}$ is

[TABLE]

where we have assumed a perturbation that preserves the rotational transform ($\delta\Phi_{1}=0$), for which $\delta B^{2}=-2\left(B^{2}\nabla\cdot\bm{\xi}_{1}+\bm{\xi}_{1}\cdot\nabla\left(B^{2}+4\pi p\right)\right)$ is the perturbation to the magnetic field strength at fixed position. The first term in (32) corresponds with the change in $f_{B}$ due to the perturbation to the field strength, while the second term accounts for the motion of the boundary. Applying the divergence theorem we obtain,

[TABLE]

The differential change in $f_{\beta}$ associated with displacement $\bm{\xi}_{1}$ satisfies

[TABLE]

The first term on the right of 34 is already in the form of a shape gradient. To evaluate the second term, we turn to the adjoint problem, and we choose $\delta\textbf{F}_{2}=-\nabla\left(\Delta p\right)$, where $\Delta\ll 1$ is a constant scalar. That is, we add a force which is proportional to the equilibrium pressure force with a small multiplier, $\Delta$. This additional force produces a proportional change in magnetic field at the boundary and thus from 22, we find

[TABLE]

Thus, we can obtain the shape gradient without perturbing the shape of the surface,

[TABLE]

Obtaining $S$ amounts to computing an equilibrium with unperturbed $\iota$, unperturbed boundary, and perturbed pressure $p^{\prime}=(1+\Delta)p$.

A similar expression can be obtained for equilibria for which the rotational transform is allowed to vary, but the toroidal current is held fixed ($\delta I_{T,1}=0$). In this case, the toroidal current for the adjoint problem is chosen to be $\delta I_{T,2}=-\Delta I_{T}\left(1/f_{P}-1/f_{B}\right)^{-1}\left(1/f_{B}\right)$ where $I_{T}$ is the unperturbed current profile, and again $\delta\textbf{F}_{2}=-\nabla\left(\Delta p\right)$. The shape gradient can then be obtained from (36).

To demonstrate, we use the NCSX stellarator LI383 equilibrium (Zarnstorff et al., 2001). The pressure profile was perturbed with $\Delta=0.01$. The unperturbed and adjoint equilibria are computed with the VMEC code (Hirshman & Whitson, 1983). The shape gradient obtained with the adjoint solution, $S_{\text{adjoint}}$, and that obtained with the direct approach, $S_{\text{direct}}$, are shown in figure 1(a). Positive values of the shape gradient indicate that $f_{\beta}$ increases if a normal perturbation is applied at a given location as indicated by (1). For the direct approach, parameter derivatives $(\partial f_{\beta}/\partial R_{mn}^{c},\partial f_{\beta}/\partial Z_{mn}^{s})$ are computed with a centered 4-point stencil for $m\leq 15$ and $|n|\leq 9$ using a polynomial fitting technique. The shape gradient is obtained using the method outlined in section 4.2 of Landreman & Paul (2018). The fractional difference between the two methods,

[TABLE]

is shown in figure 1(c). The surface-averaged value of $S_{\text{residual}}$ is $1.7\times 10^{-3}$.

The parameter $\Delta$ must be chosen carefully, as the perturbation must be large enough that the result is not dominated by round-off error, but small enough that non-linear effects do not become important. The relationship between $S_{\text{residual}}$ and $\Delta$ is shown in figure 1(d). Here $S_{\text{direct}}$ is computed using the parameters reported above such that convergence is obtained. We find that $S_{\text{residual}}$ decreases as $\left(\Delta\right)^{1}$ until $\Delta\approx 0.5$, at which point round-off error begins to dominate. This scaling is to be expected, as $\delta\textbf{B}_{2}$ is computed with a forward difference derivative with step size $\Delta$.

For this and the later examples, the computational cost of transforming the parameter derivatives to the shape gradient was negligible compared to the cost of computing the parameter derivatives. The direct approach used 2357 calls to VMEC while the adjoint approach only required two. It is clear that the adjoint method yields nearly identical derivative information to the direct method but at substantially reduced computational cost.

In figure 1 we find that the shape gradient for $f_{\beta}$ is everywhere positive. This reflects the fact that the toroidal flux enclosed by $S_{P}$ is fixed. As perturbations which displace the plasma surface outward increase the surface area of a toroidal cross-section, the toroidal field must correspondingly decrease, thus increasing $f_{\beta}$. We find that the shape gradient is increased in regions of large field strength, as indicated by the second term in (36).

3.2 Surface shape gradient for rotational transform

Consider a figure of merit, the average rotational transform in a radially localized region,

[TABLE]

Here $w(\psi)$ is a normalized weighting function,

[TABLE]

and $\psi_{m}$ and $\psi_{w}$ are parameters defining the center and width of the Gaussian weighting, respectively.

The differential change of $f_{\iota}$ associated with perturbation $\bm{\xi}_{1}$ is

[TABLE]

For the adjoint problem, we can choose a toroidal current profile $\delta I_{T,2}=I_{\Delta}w(\psi)$, where $I_{\Delta}$ is a scalar constant, and we can take $\delta\textbf{F}_{2}=0$. This additional current produces a proportional change in the magnetic field at the boundary; thus using 22, we obtain the following

[TABLE]

So, we can obtain the shape gradient from the adjoint solution

[TABLE]

Note that the computation of the shape derivative of the rotational transform on a single surface, $\psi_{m}$, with the adjoint approach would require a delta-function current perturbation, $\delta I_{T,2}=I_{\Delta}\delta(\psi-\psi_{m})$. As this type of perturbation is difficult to resolve in a numerical computation, the use of the Gaussian envelope allows the shape derivative of the rotational transform in a localized region of $\psi_{m}$ to be computed.

To demonstrate, we use the NCSX stellarator LI383 equilibrium. The toroidal current profile was perturbed with $I_{\Delta}=715$ A, $\psi_{m}=0.1\psi_{0}$, and $\psi_{w}=0.05\psi_{0}$. The shape gradient obtained with the adjoint solution and with the direct approach are shown in figure 2(a). For the direct approach, the shape gradient is computed from parameter derivatives ($\partial f_{\iota}/\partial R_{mn}^{c}$ and $\partial f_{\iota}/\partial Z_{mn}^{s}$) using a 8-point stencil with $m\leq 18$ and $|n|\leq 12$. The fractional difference, $S_{\text{residual}}$, between the two approaches is shown in figure 2(c), with a surface-averaged value of $2.7\times 10^{-2}$.

The direct approach used 7401 calls to VMEC, while the adjoint only required two. Again, it is apparent that the adjoint method allows the same derivative information to be computed at much lower computational cost.

We find that over much of the surface, the shape gradient is close to zero. A region of large negative shape gradient occurs in the concave region of the plasma surface with adjacent regions of large positive shape gradient. This indicates that “pinching” the surface in this region, making it more concave, would increase $\iota$ near the axis. The relationship between shaping and rotational transform is generally quite complex. Further analysis is needed to interpret the shape gradient for $f_{\iota}$.

3.3 Coil shape gradient for rotational transform

The shape gradient of $f_{\iota}$ can also be computed with a free boundary approach. We can again choose a toroidal current profile $\delta I_{T,2}=I_{\Delta}w(\psi)$ for the adjoint problem, where $w(\psi)$ is given by (39). Using 40 and 20 and noting that $\delta I_{T,2}$ vanishes at the boundary, we find

[TABLE]

Using 26, this can be written in terms of changes in the positions of coils in the vacuum region,

[TABLE]

When computing the coil shape gradient, the current in each coil is fixed. In arriving at 44, we assume that $\delta I_{C_{1,k}}=0$. The coil shape gradient is thus

[TABLE]

As anticipated, $\textbf{S}_{k}$ has no component in the direction tangent to the coil. Evaluating the shape gradient requires computing the adjoint magnetic field at the unperturbed coil locations in the vacuum region. This can be performed with the DIAGNO code (Gardner, 1990; Lazerson, 2012), which employs the virtual casing principle to efficiently compute the fields in the vacuum region due to the plasma current.

To demonstrate, we use the NCSX stellarator LI383 equilibrium. The toroidal current profile was perturbed with $I_{\Delta}=5.7\times 10^{5}$ A, $\psi_{m}=0.1\psi_{0}$, and $\psi_{w}=0.05\psi_{0}$. The shape gradient is computed for each of the three unique modular coils per half period (Williamson et al., 2005), keeping the planar coils fixed. The result obtained with the adjoint solution, $\textbf{S}_{\text{adjoint},k}$, is shown in figure 3. The shape gradient is also computed with the direct approach, $\textbf{S}_{\text{direct},k}$. For the direct approach, the Cartesian components of each coil are Fourier discretized ($X_{m}$, $Y_{m}$, $Z_{m})$. The numerical derivative with respect to these parameters are computed for $m\leq 45$ using an 8-point stencil. In figure 4(a) the Cartesian components of the shape gradient computed with the adjoint approach, $S_{\text{adjoint},k}^{l}$, and with the direct approach, $S_{\text{direct},k}^{l}$, are shown for each coil. Here $l\in\{x,y,z\}$. The arrows indicate the direction and magnitude of $\textbf{S}_{k}$ such that if a coil were deformed in the direction of $\textbf{S}_{k}$, $f_{\iota}$ would increase according to (2). The direct approach used 6553 calls to VMEC, while the adjoint only required two. In figure 4(b) the fractional difference between the results obtained with the two methods,

[TABLE]

is plotted. The line-averaged values of $S_{\text{residual}}^{l}$ are $6.1\times 10^{-2}$ for coil 1, $3.8\times 10^{-2}$ for coil 2, and $4.8\times 10^{-2}$ for coil 3.

From figure 3 we see that the sensitivity of $f_{\iota}$ to coil displacements is much higher in regions where the coils are close to the plasma surface. The shape gradient points toward the plasma surface in the concave region of the plasma surface, while on the outboard side the sensitivity is significantly lower, again indicating the “pinching” effect seen in figure 2.

4 Conclusions

We have obtained a relationship between 3D perturbations of MHD equilibria that is a consequence of the self-adjoint property of the MHD force operator. The relation allows for the efficient computation of shape gradients for either the outer plasma surface using the fixed boundary adjoint relation (22) or for coil shapes using the free boundary adjoint relation (20). The computation of the shape gradient of several stellarator figures of merit has been demonstrated with both the adjoint and direct approach. The application of the adjoint relation provides an $\mathcal{O}(N)$ reduction in CPU hours required in comparison with the direct method of computing the shape gradient, where $N$ is the number of parameters used to describe the shape of the outer boundary or the coils. For fully 3D geometry, $N$ can be $10^{2}-10^{3}$. Thus, the application of adjoint methods can significantly reduce the cost of computing the shape gradient for gradient-based optimization or local sensitivity analysis.

To compute the adjoint equilibria in this work, the full non-linear MHD equilibrium equations are solved using VMEC with the addition of a small perturbation to the current profile, characterized by $I_{\Delta}$, or a small perturbation to the pressure profile, characterized by $\Delta$. These parameters must be tuned carefully such that the perturbation is large enough that the result is not dominated by round-off error, but small enough that non-linear effects do not become important. If a perturbed equilibrium code were instead used, these difficulties could be avoided. However, it is convenient to use the same code for both the unperturbed and adjoint equilibrium.

It should be noted that the adjoint approach we have outlined can not yield an exact analytic shape gradient. Throughout we have assumed the existence of magnetic surfaces as the 3D equilibrium is perturbed. Therefore a code such as VMEC, which minimizes an energy subject to the constraint that surfaces exist, is suitable. Generally VMEC solutions do not satisfy (5) exactly, (Nührenberg et al., 2009), as they do not account for the formation of islands or current singularities associated with rational surfaces. Furthermore, the the parameters $\Delta$ and $I_{\Delta}$ introduce additional numerical noise. We have demonstrated that the typical difference between the shape gradient obtained with the adjoint method and that computed directly from numerical derivatives is $\lesssim 5\%$. These errors should not be significant for applying the shape gradient to an analysis of engineering tolerances. The discrepancy between the true shape gradient and that obtained numerically, with the adjoint approach or with finite difference derivatives, may become problematic as one nears a local minimum during gradient-based optimization, as the resulting shape gradient may not provide a true descent direction.

For the figure of merit considered in section 3.1, the additional force ($\delta\textbf{F}_{2}$ in (16)) applied to the adjoint problem can be expressed as a gradient of a scalar pressure; thus an existing equilibrium code could be utilized without modification. This approach could be applied to other quantities of interest. For example, consider a figure of merit which quantifies the vacuum magnetic well,

[TABLE]

where $w(\psi)$ is a weighting function and $V^{\prime\prime}(\psi)=\frac{\partial}{\partial\psi}\left(\int_{0}^{2\pi}d\theta\int_{0}^{2\pi}d\zeta\,\sqrt{g}\right)$. The presence of a magnetic well ($V^{\prime\prime}(\psi)<0$) has a stabilizing effect on MHD modes (Helander, 2014; Greene, 1997) and has been considered in stellarator design (Hirshman et al., 1999; Drevlak et al., 2018). Computing the shape gradient of $f_{W}$ requires solving an adjoint force balance equation with a perturbation to the scalar pressure, similar to the calculation in section 3.1.

For many interesting figures of merit the spatial dependence of the required force is more complicated. Thus, an equilibrium code that allows for an arbitrary force perturbation is needed. One possibility is a generalization of the code ANIMEC (Cooper et al., 2009) that currently treats anisotropic pressure tensors in the form of a bi-Maxwellian distribution. For example, the shape gradient of the neoclassical particle flux in the $1/\nu$ regime (Nemov et al., 1999) can be computed with the addition of a bulk force that takes the form $\delta\textbf{F}_{2}=-\nabla\cdot\mathbb{P}_{2}$ to the adjoint force balance equation. Here $\mathbb{P}_{2}$ is a pressure tensor that has arbitrary spatial dependence. For several figures of merit, the required additional bulk force does not take the form of the divergence of a pressure tensor. Consider the following figure of merit, which quantifies the departure from quasi-symmetry,

[TABLE]

where $w(\psi)$ is a weighting function and

[TABLE]

Here $M$ and $N$ are the mode numbers of the desired quasi-symmetry such that if $B(\psi,M\theta-N\zeta)$ for Boozer angles $\theta$ and $\zeta$, then $f_{QS}=0$ (Helander, 2014). The Boozer covariant components are $G(\psi)=2I_{P}/c$ and $I(\psi)=2I_{T}/c$, where $I_{P}$ is the poloidal current outside the $\psi$ surface. Computing the shape gradient of $f_{QS}$ requires the addition of a bulk force to the adjoint force balance equation which does not take the form of the divergence of a tensor pressure. These calculations will be reported in a separate publication. We anticipate there will be numerous additional applications of this technique for efficient optimization of MHD equilibria.

Acknowledgements

This work was supported by the US Department of Energy through grants DE-FG02-93ER-54197 and DE-FC02-08ER-54964. The computations presented in this paper have used resources at the National Energy Research Scientific Computing Center (NERSC).

Appendix A Derivation of adjoint relation

The quantity $U_{P}=U_{P_{1}}+U_{P_{2}}$ consists of two terms, accounting for changes to the vector potential due to MHD perturbations

[TABLE]

and changes to the rotational transform,

[TABLE]

The quantity $U_{P_{1}}$ can be expressed by using (16) and applying the divergence theorem to the pressure gradient terms,

[TABLE]

We will define $\delta\hat{\textbf{B}}_{1,2}=\nabla\times\left(\bm{\xi}_{1,2}\times\textbf{B}\right)$ such that $\delta\textbf{B}_{1,2}=\delta\hat{\textbf{B}}_{1,2}-\nabla\delta\Phi_{1,2}\times\nabla\zeta$. The terms in (52) due to $\delta\hat{\textbf{B}}_{1,2}$ can be evaluated using $\textbf{J}=J_{||}\textbf{b}+c\textbf{b}\times\nabla p/B$ and (5),

[TABLE]

The first term in (53) can be simplified using $\nabla\cdot\textbf{J}=0$ and noting that the perturbation can be written as $\bm{\xi}_{1,2}=\xi^{\psi}_{1,2}\nabla\psi+\xi^{\perp}_{1,2}\textbf{b}\times\nabla\psi$. Applying the identity $\textbf{b}\cdot\delta\hat{\textbf{B}}_{1,2}=-B^{2}\nabla\cdot\bm{\xi}_{1,2}-\bm{\xi}_{1,2}\cdot\nabla B^{2}-4\pi\bm{\xi}_{1,2}\cdot\nabla p$ to the second term, the following expression can be obtained,

[TABLE]

Hence we obtain the following expression for $U_{P_{1}}$,

[TABLE]

We now consider $U_{P_{2}}$ defined in (51). Applying (12) for the change in toroidal current, integrating by parts in $\psi$, and combining the expressions for $U_{P_{1}}$ (52) and $U_{P_{2}}$ (51), we obtain

[TABLE]

Next we combine $U_{P}$ (56) with $U_{B}$ (18) and add $U_{c}$ (19) to obtain the free boundary adjoint relation (20).

To obtain the fixed-boundary adjoint relation, the integral over the plasma volume (15) can be related to a surface integral by applying the divergence theorem, 21. Using (10) and applying several vector identities,

[TABLE]

Using the expression for $U_{P}$, (56), expressing the second term in (57) as a perturbed current using (12), the fixed boundary adjoint relation (22) is obtained.

Appendix B Alternate derivation of fixed-boundary adjoint relation

The MHD force operator,

[TABLE]

possesses the following self-adjointness property (Bernstein et al., 1958; Goedbloed & Poedts, 2004),

[TABLE]

where $\delta\hat{\textbf{B}}_{1,2}=\nabla\times\left(\bm{\xi}_{1,2}\times\textbf{B}\right)$ is the perturbed field corresponding to the MHD perturbations. For perturbations described by 10, 11, 12 and 16, we have the following force operator,

[TABLE]

Using 60 and several vector identities, the left hand side of 59 can be written as

[TABLE]

In arriving at (61), we use $\textbf{J}\cdot\nabla\psi=0$. Using 12 to re-express the first two terms on the right hand side,

[TABLE]

Using 10 and 59 we obtain 22.

Appendix C Interpretation of the displacement vector

For MHD perturbations such that $\delta\textbf{B}=\nabla\times\left(\bm{\xi}\times\textbf{B}\right)$ the displacement can be interpreted as a vector describing the motion of a field lines. Thus a normal perturbation to the surface of the plasma as in (1) can be expressed in terms of the displacement vector,

[TABLE]

For perturbations that allow for changes in the rotational transform it remains to be shown that a similar relation can be found.

As we require that $\psi$ remain a flux surface label in the perturbed equilibrium, the Lagrangian perturbation to $\psi$ at fixed position is

[TABLE]

The perturbed magnetic field, $\textbf{B}^{\prime}=\textbf{B}+\delta\textbf{B}$ must remain tangent to $\psi^{\prime}=\psi+\delta\psi$ surfaces; thus to first order in the perturbation,

[TABLE]

Applying the form for the perturbed field allowing for changes in the rotational transform, $\delta\textbf{B}=\nabla\times\left(\bm{\xi}\times\textbf{B}-\delta\Phi\nabla\zeta\right)$, and using several vector identities, the following condition is obtained

[TABLE]

This implies that $\delta\textbf{r}\cdot\nabla\psi=\bm{\xi}\cdot\nabla\psi+F(\psi)$, where $F(\psi)$ is some flux function which can be determined by requiring that the perturbation to the toroidal flux as a function of $\psi$ vanishes, $\delta\Phi_{T}(\psi)=0$.

The perturbed toroidal flux through a surface labeled by $\psi$ contains two terms, corresponding to the flux of the unperturbed field through the perturbed surface and the perturbed field through the unperturbed surface,

[TABLE]

Using the form for $\delta\textbf{B}$, applying the divergence theorem, and noting that $\textbf{B}\cdot\nabla\zeta=\sqrt{g}^{-1}$, the following condition is obtained,

[TABLE]

By requiring that $\delta\Phi_{T}(\psi)=0$, we find that $F(\psi)=0$. Thus we can express shape gradients in the form of (63) even when the rotational transform is allowed to vary.