2. Basics of adjoint methods

This section aims to briefly introduce adjoint methods, without delving deeply in the mathematical formalism. More details can be found in, for example, Delfour & Zolésio (Reference Delfour and Zolésio2011) and Paul et al. (Reference Paul, Antonsen, Landreman and Cooper2020). Note that here we exclusively consider volumes $\mathcal {V}$ with closed, smooth, connected and orientable boundaries $\mathcal {S}=\partial \mathcal {V}$, such that there is a one-to-one correspondence between $\mathcal {V}$ and $\mathcal {S}$, a proof of which can be found in McGrath (Reference McGrath2016). For the specific case of closed $\mathcal {S}$, it is therefore equivalent to denote shape functionals as depending either on the volume, e.g. $J(\mathcal {V})$, or on the surface, e.g. $J(\mathcal {S})$. We choose the latter approach for greater clarity to the reader, as the boundary is the object that will ultimately be optimised over. However, note that the mathematical literature (Delfour & Zolésio Reference Delfour and Zolésio2011) would emphasise the dependence on volumes instead.

We are interested in obtaining derivative information for a shape functional $J(\mathcal {S}) = f(\mathcal {S}, u(\mathcal {S}))$, called hereafter the objective function. This functional depends on the closed surface $\mathcal {S}$ explicitly and also implicitly through the solution $u(\mathcal {S})$ to a partial differential equation (PDE) $\mathcal {P}(\mathcal {S},u(\mathcal {S})) = 0$ on $\mathcal {S}$. Here, we consider $u$ to be a scalar function on the boundary $\mathcal {S}$, with a sufficient number of continuous derivatives for the operator $\mathcal {P}$ under consideration. The differential operator $\mathcal {P}$ can be linear or nonlinear, and is assumed to map back to a space of scalar functions on the boundary.

Consider a displacement of the surface $\mathcal {S}$ in the direction $\delta \boldsymbol {x}$ with magnitude $\epsilon$, resulting in a perturbed surface $\mathcal {S}_\epsilon = \{ \boldsymbol {x}_0 + \epsilon \delta \boldsymbol {x}(\boldsymbol {x}_0) : \boldsymbol {x}_0 \in \mathcal {S}\}$. The shape derivative of an arbitrary function $g(\mathcal {S})$ in the direction $\delta \boldsymbol {x}$ is now defined as

(2.1)\begin{equation} \delta g(\mathcal{S})[\delta\boldsymbol{x}] = \lim_{\epsilon\rightarrow 0} \frac{g(\mathcal{S_\epsilon}) - g(\mathcal{S})}{\epsilon}. \end{equation}

If $g(\mathcal {S})$ depends only on the geometrical shape of the surface, the shape derivative $\delta g(\mathcal {S})[\delta \boldsymbol {x}]$ will be a function of only the normal component $ {\delta \boldsymbol {{x}}\boldsymbol {\cdot } {\boldsymbol {\hat {{n}}}}}$ of the displacement, as any tangential component of $\delta \boldsymbol {x}$ leaves the shape of $\mathcal {S}$ unchanged to first order. Here, $ {\boldsymbol {\hat {{n}}}}$ is a normal vector on $\mathcal {S}$.

Obtaining the shape derivative of the objective function is non-trivial due to the implicit dependence on $\mathcal {S}$ through $u(\mathcal {S})$. Let us consider instead the Lagrangian

(2.2)\begin{equation} \mathcal{L}(\mathcal{S},u,q) = f(\mathcal{S},u) + \int_\mathcal{S} \,\mathrm{d} S \; q \; \mathcal{P}(\mathcal{S},u), \end{equation}

with the Lagrange multiplier $q$. Note that $\mathcal {S}$, $u$ and $q$ are independent variables in the above definition. When $u=u(\mathcal {S})$ such that the PDE constraint $\mathcal {P}(\mathcal {S},u(\mathcal {S})) = 0$ is satisfied, the Lagrangian is equal to the objective function for any choice of $q$, i.e. $\mathcal {L}(\mathcal {S}, u(\mathcal {S}), q)=f(\mathcal {S}, u(\mathcal {S})) = J(\mathcal {S})$. It then follows from (2.1) that $\delta \mathcal {L}(\mathcal {S}, u(\mathcal {S}), q)[\delta \boldsymbol {x}]=\delta J(\mathcal {S})[\delta \boldsymbol {x}]$, i.e. we can obtain the objective function's shape derivative by computing that of the Lagrangian. The adjoint method now makes use of the freedom in $q$ to cancel the contribution $\delta u(\mathcal {S})[\delta \boldsymbol {x}]$ contained in $\delta \mathcal {L}(\mathcal {S}, u(\mathcal {S}), q)[\delta \boldsymbol {x}]$, i.e. $q$ is chosen such that the Lagrangian is stationary with respect to $u(\mathcal {S})$. This leads to an adjoint PDE for $q$, the solution of which we denote by $q(\mathcal {S})$.

Due to the Hadamard–Zolésio structure theorem (Delfour & Zolésio Reference Delfour and Zolésio2011), the shape derivative of our objective function can ultimately be expressed as

(2.3)\begin{equation} \delta J(\mathcal{S})[\delta\boldsymbol{x}] = \delta\mathcal{L}(\mathcal{S},u(\mathcal{S}),q(\mathcal{S}))[\delta\boldsymbol{x}] = \int_\mathcal{S} \,\mathrm{d} S \;({\delta\boldsymbol{{x}}\boldsymbol{\cdot} {\boldsymbol{\hat{{n}}}}})\; \mathcal{G}(\mathcal{S},u(\mathcal{S}),q(\mathcal{S})), \end{equation}

where $\mathcal {G}$ is called the shape gradient, and can be interpreted as the local sensitivity of the objective function to perturbations of $\mathcal {S}$.

In practice, the surface $\mathcal {S}$ is typically represented by a finite set of parameters $\varOmega = \{\varOmega _i,\; i = 1, 2, \ldots N\}$, e.g. Fourier coefficients $\{R_{m,n}, Z_{m,n}\}$. The functional $J(\mathcal {S})$ is approximated by a function $J(\varOmega )$, and similarly for $\mathcal {G}$, $u(\mathcal {S})$ and $q(\mathcal {S})$. The derivative of $J(\varOmega )$ with respect to a parameter $\varOmega _i$ is then approximated as

(2.4)\begin{equation} \frac{\partial J(\varOmega)}{\partial \varOmega_i} = \int_\mathcal{S} \,\mathrm{d} S\; \frac{\partial \boldsymbol{x}}{\partial \varOmega_i}\boldsymbol{\cdot} {\boldsymbol{\hat{{n}}}} \; \mathcal{G}(\varOmega,u(\varOmega),q(\varOmega)). \end{equation}

When evaluating the parameter derivatives numerically through (2.4), errors are introduced from the inexact solutions to the original and adjoint PDEs. Indeed, these PDEs are assumed to be exactly satisfied in the preceding derivation, and also typically when deriving an expression for the shape gradient $\mathcal {G}$.

The formalism presented above can be generalised to multiple PDE constraints. We can also consider PDEs satisfied not on $\mathcal {S}$ but in the volume enclosed by it, in which case $u$ would be a scalar function in the volume, the second term on the right-hand side of (2.2) would be a volume integral and the Hadamard–Zolésio structure theorem (Delfour & Zolésio Reference Delfour and Zolésio2011) would still guarantee that we can ultimately express the shape derivative of the objective as the boundary integral (2.3). In § 4, we consider two PDE constraints: the Laplace equation for the vacuum field and the straight field line equation, respectively valid in the volume and on the boundary, with two corresponding adjoint variables.

3. Evaluating approximate flux coordinates on an isolated flux surface

The existence of toroidally nested flux surfaces (closed orientable surfaces on which the magnetic field is everywhere tangent and non-vanishing) is commonly assumed in theoretical studies of magnetically confined plasmas. In particular, many formulas involve $\boldsymbol {\nabla }\psi$, where the toroidal flux $\psi$ is a global flux surface label. However, three-dimensional magnetic fields lacking a continuous symmetry are not generally integrable. It is desirable to generalise $\boldsymbol {\nabla }\psi$ to the case of an isolated flux surface, i.e. a flux surface in whose neighbourhood the field is generally non-integrable.

Consider a toroidal flux surface $\mathcal {S}$ with a magnetic field that is nowhere vanishing. On $\mathcal {S}$, the magnetic field's normal component vanishes by definition, i.e. $\boldsymbol {B}\boldsymbol {\cdot } {\boldsymbol {\hat {{n}}}} = 0$ with $ {\boldsymbol {\hat {{n}}}}$ the unit normal vector on $\mathcal {S}$. The field line label $\alpha$ on $\mathcal {S}$ is defined through the straight field line equation $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }_{\varGamma }\alpha = 0$. Here, the tangential gradient $\boldsymbol {\nabla }_{\varGamma }$, defined in Appendix A.1, is the component of the three-dimensional gradient tangential to the surface (A 1). Uniqueness of $\boldsymbol {\nabla }_{\varGamma }\alpha$ is guaranteed by fixing the secular component in the poloidal angle $\theta$ to be unity, i.e. $\alpha =\theta -\iota \phi +\lambda (\theta,\phi )$, with arbitrary poloidal and toroidal angles $\theta$ and $\phi$ on $\mathcal {S}$, the rotational transform $\iota$ and the single-valued function $\lambda$. We here assume the toroidal component of the field $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\phi$ is nowhere vanishing, as this allows the construction of an area-preserving mapping, from a toroidal cut of the boundary onto itself (Meiss Reference Meiss1992). This mapping can be associated with a rotation number (Greene, Mackay & Stark Reference Greene, Mackay and Stark1986) equivalent to the rotational transform $\iota$, implying that the straight field line equation can be solved.

In the integrable case, $\boldsymbol {\nabla }\psi$ is normal to the flux surfaces, and the magnetic field satisfies $\boldsymbol {B} = \boldsymbol {\nabla }\psi \times \boldsymbol {\nabla }\alpha$. We now define the generalisation $\boldsymbol {g}_{\psi }$ on $\mathcal {S}$ of the toroidal flux gradient $\boldsymbol {\nabla }\psi$, by setting $\boldsymbol {g}_{\psi }$ normal to $\mathcal {S}$, and by requiring $\boldsymbol {B} = \boldsymbol {g}_{\psi }\times \boldsymbol {\nabla }\alpha$ to be satisfied on $\mathcal {S}$. Squaring the latter equality and using $\boldsymbol {g}_{\psi } = {\boldsymbol {\hat {{n}}}} |\boldsymbol {g}_{\psi }|$ yields

(3.1)\begin{equation} \boldsymbol{g}_{\psi} = {\boldsymbol{\hat{{n}}}} \; \frac{B}{{\left| \boldsymbol{\nabla}_{\varGamma} \alpha \right|}}, \end{equation}

where $B = {\left | \boldsymbol {B} \right |}$ is the magnetic field strength.

The defining expression for $\boldsymbol {g}_{\psi }$ (3.1) can be evaluated on any flux surface without requiring nested flux surfaces in its neighbourhood, and will revert to $\boldsymbol {g}_{\psi } = \boldsymbol {\nabla }\psi$ when the field is integrable in the neighbourhood of that flux surface. In practice, one might couple objective functions relying on (3.1) with a figure of merit targeting integrability, aiming for a final plasma shape for which the field is integrable, such that $\boldsymbol {g}_{\psi } = \boldsymbol {\nabla }\psi$ and the minimised objective function represents the physical quantity of interest.

The generalised toroidal flux gradient (3.1) is evaluated in figure 1(b) for a rotating ellipse configuration, with the magnetic field computed using the SPEC code and the straight field line equation solved using a Fourier–Galerkin spectral solver. It agrees excellently with the toroidal flux gradient evaluated by VMEC, which can be calculated directly due to the imposed nestedness of flux surfaces, shown in figure 1(a). The relative difference between the two is below a percentage point in this case, as shown in figure 1(c). The difference is expected to be small when integrability is satisfied, which indeed seems to hold here, as attested by the absence of islands and chaotic regions in the SPEC Poincaré plot shown in figure 1(d). Note that SPEC solves for a vacuum magnetic field, while VMEC computes an ideal MHS equilibrium with vanishing thermal pressure, and with a plasma current that is small but finite due to the constraint of integrability.

Figure 1. Comparison of (a) the toroidal flux gradient ${\left | \boldsymbol {\nabla }\psi \right |}$ evaluated using VMEC with (b) the generalised toroidal flux gradient $|\boldsymbol {g}_{\psi }|$ (3.1) obtained in SPEC, for a $5$-period rotating ellipse case with half a rotation per field period, major radius at the ellipse centre $R_0 = 5$ m and ellipse major and minor axes values of $2$ and $1$ m, respectively. The relative difference between the two quantities is below $1\,\%$, as shown in (c). A small difference is to be expected in this case, where integrability is well satisfied, as attested in (d) by the Poincaré plot at toroidal angle $\phi =0$ from the SPEC calculation, which agrees well with the flux surfaces computed by VMEC. All data generated in this paper can be obtained from Nies (Reference Nies2021).

The generalised toroidal flux gradient $\boldsymbol {g}_{\psi }$ can be applied generally in any situation where local flux coordinates need to be evaluated, e.g. in calculations of perpendicular transport or magnetohydrodynamic stability. Isolated flux surfaces occur, for example, in fixed-boundary equilibrium calculations, where the plasma outer boundary is constrained to be a flux surface as a boundary condition on the magnetic field, or at the interfaces of multi-region relaxed magnetohydrodynamic (MRxMHD) equilibria computed by, for example, SPEC (Hudson et al. Reference Hudson, Dewar, Dennis, Hole, McGann, von Nessi and Lazerson2012) or BIEST (Malhotra et al. Reference Malhotra, Cerfon, Imbert-Gérard and O'Neil2019). In the following, we employ (3.1) specifically for a fixed-boundary vacuum field to formulate quasi-symmetry on the boundary.

4. Application of adjoint formalism to vacuum fields

Consider a vacuum magnetic field $\boldsymbol {B}$ in a toroidal domain $\mathcal {V}$ bounded by the surface $\mathcal {S} = \partial \mathcal {V}$. As the vacuum magnetic field is curl-free, it can be expressed as $\boldsymbol {B}=\boldsymbol {\nabla }\varPhi$, with the scalar potential $\varPhi$. Because we consider a simple torus $\mathcal {V}$, the most general form for the scalar potential is $\varPhi =G (\omega +\phi )$, where $G$ is a constant, $\omega$ is a single-valued function on $\mathcal {S}$ and $\phi$ is an arbitrary toroidal angle. By integrating the magnetic field along a toroidal loop around the torus, the constant $G$ is found to be proportional to the net external current through the ‘hole’ of the torus.

As the magnetic field is divergence-less, the magnetic scalar potential satisfies the Laplace equation. The field's normal component is constrained to vanish on $\mathcal {S}$ by imposing a Neumann boundary condition on the magnetic scalar potential. Further prescribing, for example, $G$, or the toroidal flux, guarantees a unique solution to Laplace's equation. We herein opt to hold the toroidal flux fixed, although the shape derivative $\delta G[\delta \boldsymbol {x}]$ will not appear in this study due to our normalisation of the figure of merit for quasi-symmetry (4.18). A different choice of normalisation would lead to an additional contribution proportional to $\delta G[\delta \boldsymbol {x}]$ in the shape derivative of the Lagrangian.

For convenience, we define the normalised magnetic field $\boldsymbol {\breve {B}}$ as

(4.1)\begin{equation} \boldsymbol{\breve{B}} \equiv \boldsymbol{B}/G = \boldsymbol{\nabla}\left(\omega + \phi\right). \end{equation}

Let us further assume the toroidal angle $\phi$ to be the azimuthal angle in cylindrical coordinates, satisfying $\Delta \phi =0$ in the domain of interest. We can thus write

(4.2a)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\breve{B}} = \Delta\omega = 0, \quad \mathrm{in}\ \mathcal{V}, \end{gather}

(4.2b)\begin{gather}\boldsymbol{\breve{B}}\boldsymbol{\cdot}{\boldsymbol{\hat{{n}}}} = \boldsymbol{\nabla}(\omega+\phi)\boldsymbol{\cdot}{\boldsymbol{\hat{{n}}}} = 0, \quad\mathrm{on}\ \mathcal{S}, \end{gather}

with $ {\boldsymbol {\hat {{n}}}}$ the normal unit vector on $\mathcal {S}$. Furthermore, the shape derivative $\delta \omega [\delta \boldsymbol {x}]$ satisfies

(4.3a)\begin{gather} \Delta(\delta\omega[\delta\boldsymbol{x}]) = 0, \quad \mathrm{in}\ \mathcal{V}, \end{gather}

(4.3b)\begin{gather}\boldsymbol{\breve{B}}\boldsymbol{\cdot}\delta{\boldsymbol{\hat{{n}}}}[\delta\boldsymbol{x}] + \boldsymbol{\nabla}(\delta\omega[\delta\boldsymbol{x}])\boldsymbol{\cdot}{\boldsymbol{\hat{{n}}}} + ({\delta\boldsymbol{{x}}\boldsymbol{\cdot} {\boldsymbol{\hat{{n}}}}}){\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{\breve{B}}\boldsymbol{\cdot}{\boldsymbol{\hat{{n}}}} = 0, \quad\mathrm{on}\ \mathcal{S}, \end{gather}

where the Laplace equation is obtained from noting the commutative property of shape and spatial derivatives, and the normal boundary condition on $\delta \omega$ was derived in, for example, Sokołowski & Zolésio (Reference Sokołowski and Zolésio1992, § 3.2). The shape derivative of the normal vector $\delta {\boldsymbol {\hat {{n}}}}[\delta \boldsymbol {x}]=-\boldsymbol {\nabla }_{\varGamma }( {\delta \boldsymbol {{x}}\boldsymbol {\cdot } {\boldsymbol {\hat {{n}}}}})$ is derived in Appendix A.3.

Evaluating the rotational transform and quasi-symmetry figures of merit further requires the solution to the straight field line equation

(4.4)\begin{equation} 0 = \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\varGamma}\alpha, \quad\mathrm{on}\ \mathcal{S}, \end{equation}

with the field line label

(4.5)\begin{equation} \alpha \equiv \theta-\iota\phi+\lambda(\theta,\phi), \end{equation}

where $\theta$ is a general poloidal angle, $\lambda$ is a single-valued function of $\theta$ and $\phi$ and $\iota$ is a scalar. Note that both $\lambda$ and $\iota$ are defined on the boundary $\mathcal {S}$ only, through (4.5).

Let us define the Lagrangian corresponding to an arbitrary objective function $f(\mathcal {S},\omega,\iota,\lambda )$:

(4.6)\begin{equation} \mathcal{L}(\mathcal{S}, \omega, q_\omega, \iota, \lambda, q_\alpha ) = f(\mathcal{S},\omega,\iota,\lambda) + \mathcal{M}(\mathcal{S}, \omega, q_\omega) + \mathcal{N}(\mathcal{S}, \omega, \iota, \lambda, q_\alpha), \end{equation}

with the weak form of the Laplace equation (4.2a):

(4.7)\begin{equation} \mathcal{M}(\mathcal{S}, \omega, q_\omega) = \int_{\mathcal{V}} \,\mathrm{d} V \; q_\omega \Delta \omega \end{equation}

and the weak form of the straight field line equation (4.4) normalised by $G$:

(4.8)\begin{equation} \mathcal{N}(\mathcal{S}, \omega, \iota, \lambda, q_\alpha) = \int_{\mathcal{S}} \,\mathrm{d} S \; q_\alpha \boldsymbol{\breve{B}} \boldsymbol{\cdot} \boldsymbol{\nabla}_{\varGamma} \alpha. \end{equation}

Here, the Lagrange multipliers $q_\alpha$ and $q_\omega$ are single-valued functions on the surface $\mathcal {S}$ and in the volume $\mathcal {V}$, respectively.

In the following, we assume that $\omega =\omega (\mathcal {S})$ satisfies (4.2a) and $\alpha =\alpha (\mathcal {S})$ satisfies (4.4). For ease of notation, we do not write out the dependencies on $\mathcal {S}$ explicitly, and we also omit writing the arguments of, for example, $f$, $\mathcal {M}$, $\mathcal {N}$ and $\mathcal {L}$. This means that the shape derivative $\delta \mathcal {M}[\delta \boldsymbol {x}]$ (4.9) should be understood as $\delta \mathcal {M}(\mathcal {S}, \omega (\mathcal {S}), q_\omega )[\delta \boldsymbol {x}]$, with contributions from $\delta \omega (\mathcal {S})[\delta \boldsymbol {x}]$; and similarly for other shape derivatives.

To evaluate the shape derivative of the Lagrangian (4.6), we first compute the shape derivatives of $\mathcal {M}$ and $\mathcal {N}$, before turning to the shape derivatives of our objective functions in §§ 4.1 and 4.2. First, the shape derivative of $\mathcal {M}$ (4.7), derived in Appendix B.1, is

(4.9)\begin{equation} \delta\mathcal{M}[\delta\boldsymbol{x}] = \int_{\mathcal{V}} \,\mathrm{d} V \; \delta\omega[\delta\boldsymbol{x}] \Delta q_\omega - \int_{\mathcal{S}} \,\mathrm{d} S \; \left[ \delta \omega[\delta\boldsymbol{x}] \boldsymbol{\nabla} q_\omega \boldsymbol{\cdot} {\boldsymbol{\hat{{n}}}} - ({\delta\boldsymbol{{x}}\boldsymbol{\cdot} {\boldsymbol{\hat{{n}}}}}) \boldsymbol{\breve{B}}\boldsymbol{\cdot}\boldsymbol{\nabla} q_\omega \right]. \end{equation}

Second, the shape derivative of $\mathcal {N}$ (4.8), derived in Appendix B.2, is

(4.10)\begin{align} \delta\mathcal{N}[\delta\boldsymbol{x}] & = \int_{\mathcal{S}} \mathrm{d} S \left[ - \delta\omega[\delta\boldsymbol{x}] \boldsymbol{\nabla}_{\varGamma} \boldsymbol{\cdot} \left(q_\alpha \boldsymbol{\nabla}_{\varGamma} \alpha\right) - \delta\iota[\delta\boldsymbol{x}] q_\alpha \boldsymbol{\breve{B}} \boldsymbol{\cdot} \boldsymbol{\nabla}\phi\right.\nonumber\\ & \quad \left.- \delta\lambda[\delta\boldsymbol{x}] \boldsymbol{\nabla}_{\varGamma} \boldsymbol{\cdot} \left( q_\alpha \boldsymbol{\breve{B}} \right) + ({\delta\boldsymbol{{x}}\boldsymbol{\cdot} {\boldsymbol{\hat{{n}}}}}) q_\alpha \left( {\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\breve{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\varGamma}\alpha - \boldsymbol{\breve{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\varGamma}\alpha \right) \right]. \end{align}

The tangential gradient $\boldsymbol {\nabla }_{\varGamma }(\boldsymbol {\cdot })$ and tangential divergence $\boldsymbol {\nabla }_{\varGamma }\boldsymbol {\cdot }(\boldsymbol {\cdot })$ operators are defined in Appendix A.1.

We now proceed by computing the shape derivatives of two objective functions, first targeting a given rotational transform value on $\mathcal {S}$ (§ 4.1) and second targeting quasi-symmetry on $\mathcal {S}$ with a given helicity value (§ 4.2). We will then be able to evaluate the shape derivative of the Lagrangian (4.6), yielding the adjoint equations and shape gradient formulas. Numerical verification and example shape gradients are shown for each figure of merit.

4.1. Rotational transform objective function

Before evaluating the more complicated shape gradient for the quasi-symmetry figure of merit in § 4.2, we consider a simple figure of merit targeting a given target rotational transform $\iota _T$ on the surface $\mathcal {S}$. We thus define

(4.11)\begin{equation} f_\iota = \tfrac{1}{2} ( \iota - \iota_T )^2, \end{equation}

where $\iota$ is the rotational transform on $\mathcal {S}$, obtained by solving the straight field line equation (4.4). The shape derivative of $f_\iota$ is simply

(4.12)\begin{equation} \delta f_\iota [\delta\boldsymbol{x}] = \delta \iota[\delta\boldsymbol{x}] ( \iota - \iota_T ). \end{equation}

By combining (4.9), (4.10) and (4.12), we obtain the shape derivative of the Lagrangian $\mathcal {L}_\iota$ ((4.6) with $f=f_\iota$):

(4.13)\begin{align} \delta\mathcal{L}_\iota[\delta\boldsymbol{x}] & = \int_{\mathcal{V}} \,\mathrm{d} V \; \delta\omega[\delta\boldsymbol{x}] \Delta q_\omega - \int_{\mathcal{S}} \,\mathrm{d} S \; \delta \omega [\delta\boldsymbol{x}] \left[ \boldsymbol{\nabla} q_\omega \boldsymbol{\cdot} {\boldsymbol{\hat{{n}}}} + \boldsymbol{\nabla}_{\varGamma} \boldsymbol{\cdot} \left(q_\alpha \boldsymbol{\nabla}_{\varGamma} \alpha \right) \right]\nonumber\\ & \quad - \delta\iota[\delta\boldsymbol{x}] \left[ \int_{\mathcal{S}} \,\mathrm{d} S \; q_\alpha \boldsymbol{\breve{B}} \boldsymbol{\cdot} \boldsymbol{\nabla}\phi - ( \iota - \iota_T ) \right] - \int_{\mathcal{S}} \,\mathrm{d} S \;\delta\lambda[\delta\boldsymbol{x}] \boldsymbol{\nabla}_{\varGamma} \boldsymbol{\cdot} \left( q_\alpha \boldsymbol{\breve{B}} \right)\nonumber\\ & \quad + \int_{\mathcal{S}} \,\mathrm{d} S ({\delta\boldsymbol{{x}}\boldsymbol{\cdot} {\boldsymbol{\hat{{n}}}}}) \left[\boldsymbol{\nabla} q_\omega \boldsymbol{\cdot} \boldsymbol{\breve{B}} + q_\alpha \left({\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\breve{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\varGamma}\alpha - \boldsymbol{\breve{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\varGamma}\alpha \right) \right]. \end{align}

First, we obtain the adjoint equation for $q_\alpha$ by requiring the second line of (4.13) to vanish:

(4.14a)\begin{gather} \boldsymbol{\nabla}_{\varGamma} \boldsymbol{\cdot} \left( q_\alpha \boldsymbol{\breve{B}} \right) = 0, \end{gather}

(4.14b)\begin{gather}\int_{\mathcal{S}} \,\mathrm{d} S \; q_\alpha \boldsymbol{\breve{B}} \boldsymbol{\cdot} \boldsymbol{\nabla}\phi - ( \iota - \iota_T ) = 0, \end{gather}

with both equations defined on $\mathcal {S}$. Using (A 3), the surface integral of (4.14a) yields $0 = \boldsymbol {\breve {B}}\boldsymbol {\cdot } {\boldsymbol {\hat {{n}}}}$, which is consistent with the boundary condition on the magnetic field (4.2b). The first equation (4.14a) can be recast in the form of a magnetic differential equation $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla } q_\alpha = - q_\alpha \boldsymbol {\nabla }_{\varGamma }\boldsymbol {\cdot }\boldsymbol {B}$, while the second equation (4.14b) is an integral condition on $q_\alpha$ that ensures uniqueness of the solution.

Second, the adjoint equation for $q_\omega$ is obtained by requiring the first line of (4.13) to vanish:

(4.15a)\begin{gather} \Delta q_\omega = 0, \quad \mathrm{in}\ \mathcal{V}, \end{gather}

(4.15b)\begin{gather}\boldsymbol{\nabla} q_\omega\boldsymbol{\cdot}{\boldsymbol{\hat{{n}}}} ={-}\boldsymbol{\nabla}_{\varGamma} \boldsymbol{\cdot} \left(q_\alpha \boldsymbol{\nabla}_{\varGamma} \alpha \right), \quad\mathrm{on}\ \mathcal{S}. \end{gather}

Like the magnetic potential $\omega$, the adjoint variable $q_\omega$ satisfies the Laplace equation in $\mathcal {V}$ (4.15a). However, contrary to $\omega$, $q_\omega$ has a non-zero normal boundary condition on $\mathcal {S}$ (4.15b), which notably depends on the straight field line adjoint variable $q_\alpha$. Equations (4.15b) and (4.15a) are consistent, as $\int _\mathcal {V} \,\mathrm {d} V\; \Delta q_\omega = \int _\mathcal {S} \,\mathrm {d} S\; \boldsymbol {\nabla } q_\omega \boldsymbol {\cdot } {\boldsymbol {\hat {{n}}}} = 0$, by (A 3).

Finally, the remaining contribution from the last line of (4.13) yields the shape gradient

(4.16)\begin{equation} \mathcal{G_\iota} = \frac{1}{G} \left[ \boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{\nabla} q_\omega + q_\alpha \left( {\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{B} - \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{\hat{{n}}}} \right)\boldsymbol{\cdot}\boldsymbol{\nabla}_{\varGamma} \alpha \right], \end{equation}

with $\delta \mathcal {L}_\iota [\delta \boldsymbol {x}] = \int _{\mathcal {S}} \,\mathrm {d} S ( {\delta \boldsymbol {{x}}\boldsymbol {\cdot } {\boldsymbol {\hat {{n}}}}}) \mathcal {G_\iota }$.

We now calculate the shape gradient numerically and verify it against a finite-difference evaluation. The solutions to Laplace's equation for the vacuum magnetic field (4.2a) and adjoint equation for $q_\omega$ (4.15a) are calculated with the SPEC code (Hudson et al. Reference Hudson, Dewar, Dennis, Hole, McGann, von Nessi and Lazerson2012), employing the new Zernike polynomial implementation (Qu et al. Reference Qu, Pfefferlé, Hudson, Baillod, Kumar, Dewar and Hole2020). In all results shown, the radial resolution $L_\mathrm {rad}$ in SPEC is tied to the poloidal Fourier resolution $M_\mathrm {pol}$ through $L_\mathrm {rad} = M_\mathrm {pol} + 4$. The solutions to the straight field line and $q_\alpha$ adjoint equations are obtained with a Fourier–Galerkin spectral solver.

The shape gradient (4.16) is shown in figure 2(a) for the example rotating ellipse case introduced in figure 1. The localisation at the ellipse tips is unsurprising, as near-axis expansions show that ellipticity of the flux surfaces generates rotational transform (Mercier Reference Mercier1964). This shape gradient $\mathcal {G}_\mathrm {adjoint}$ can be verified against the direct finite-difference evaluation $\mathcal {G}_\mathrm {FD}$ shown in figure 2(b), obtained by evaluating the parameter derivatives $\partial f/\partial \varOmega _i$ through finite differences and inverting (2.4) (see Landreman & Paul Reference Landreman and Paul2018). On the scale of the figure, the two shape gradients seem identical. The relative error is shown in figure 2(c) to be small, limited to ${\sim }2\,\%$ at the ellipse tips, and exhibits oscillations typical of a truncated Fourier resolution. The relative error is here defined as the absolute error normalised by the $L^\infty$-norm of $\mathcal {G}_\mathrm {adjoint}$, i.e. its maximum absolute value. This choice is preferable to, for example, the $L^2$-norm, as the shape gradient and the error thereof have small average values on the boundary compared withtheir large values at the ellipse tips, such that unreasonably high relative errors would result at these locations if using the $L^2$-norm as normalisation.

Figure 2. Shape gradient for the rotational transform objective function with $\iota _T = 1$, evaluated through (a) adjoint methods and (b) a forward finite-difference scheme with step size $\epsilon _\mathrm {FD}=10^{-7}$, for the example rotating ellipse case introduced in figure 1, with Fourier resolution $(N_\mathrm {tor}, M_\mathrm {pol}) = (16,16)$. The relative error, defined as the absolute error normalised by the maximal absolute value of the adjoint shape gradient, is shown in (c). The convergence of the relative error in the parameter derivative (2.4) for a random direction in $\varOmega$ is shown in (d) as a function of the step size $\epsilon _\mathrm {FD}$ and Fourier resolution $(N_\mathrm {tor}, M_\mathrm {pol})$. The black dashed line indicates the linear scaling in $\epsilon _\mathrm {FD}$ expected from the employed forward finite-difference scheme.

Furthermore, we test convergence of the shape gradient by evaluating a parameter derivative (2.4) for a random direction in $\varOmega$. The parameter derivative is evaluated both through the adjoint shape gradient and by a forward finite-difference scheme. The relative error is shown in figure 2(d) as a function of the finite-difference step size $\epsilon _\mathrm {FD}$ and the Fourier resolution, which is used in both SPEC and the Fourier–Galerkin spectral solver. As $\epsilon _\mathrm {FD}$ is reduced, the error initially decreases linearly with $\epsilon _\mathrm {FD}$, as expected from the employed forward finite-difference scheme, until it plateaus at a value governed by the finite Fourier or radial resolution. As mentioned in § 2, errors in the adjoint shape gradient are introduced by the assumption that the constraint and adjoint PDEs are exactly satisfied. In practice, these PDEs are solved only approximately, limited by the finite Fourier and radial resolution, such that a reduction of the error with increasing resolution is to be expected.

4.2. Quasi-symmetry objective function

For a general (non-vacuum) magnetic field with nested flux surfaces, quasi-symmetry can be expressed as

(4.17)\begin{equation} \frac{\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}\psi\times\boldsymbol{\nabla} B}{\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} B} ={-}\frac{MG + NI}{N-\iota M}, \end{equation}

where $I$ is the net toroidal plasma current and $N/M$ is the helicity of the field strength in Boozer coordinates (see e.g. Helander Reference Helander2014). For the vacuum field considered here, $I=0$. In the following, we do not consider quasi-poloidal symmetry, i.e. we assume $M\ne 0$. If desired, it would be straightforward to extend the derived results to include the case $M=0$.

For magnetic fields with globally nested flux surfaces labelled by $\psi$, (4.17) is defined globally. However, we are considering a generally non-integrable field, assuming only that the boundary $\mathcal {S}$ is a flux surface. Using the generalised toroidal flux gradient defined in (3.1), we are able to define quasi-symmetry on the isolated flux surface $\mathcal {S}$, leading to the quasi-symmetry objective function

(4.18)\begin{equation} f_\mathrm{QS} = \frac{1}{2} \int_\mathcal{S} \,\mathrm{d} S\; v_\mathrm{QS}^2, \end{equation}

with

(4.19)\begin{equation} v_\mathrm{QS} = \boldsymbol{\breve{B}}\boldsymbol{\cdot}\boldsymbol{\nabla} \breve{B} - \boldsymbol{\breve{B}}\times\boldsymbol{\breve{g}}_\psi\boldsymbol{\cdot}\boldsymbol{\nabla} \breve{B} \left( \iota - N/M\right), \end{equation}

where we defined $\boldsymbol {\breve {g}}_\psi \equiv \boldsymbol {g}_{\psi }/G$. If $f_\mathrm {QS}=0$ and the field is integrable in the neighbourhood of $\mathcal {S}$, (4.17) will be satisfied on $\mathcal {S}$, i.e. the field is quasi-symmetricon the boundary.

The shape derivative of $f_\mathrm {QS}$ is derived in Appendix B.3, with the final expression given in (B 24). Combined with the shape derivatives of $\mathcal {M}$ (4.9) and $\mathcal {N}$ (4.10), the shape derivative of the Lagrangian (4.6) with the quasi-symmetric figure of merit follows (B 26).

Requiring the Lagrangian to be stationary with respect to variations in $\iota$ and $\lambda$, the first two lines of (B 26) yield the adjoint equations for $q_\alpha$:

(4.20a)\begin{gather} \boldsymbol{\nabla}_{\varGamma} \boldsymbol{\cdot} \left( q_\alpha \boldsymbol{\breve{B}} \right) ={-}\boldsymbol{\nabla}_{\varGamma}\boldsymbol{\cdot}\left[ \boldsymbol{\nabla}_{\varGamma}\alpha \left( v_\mathrm{QS}\; \boldsymbol{\breve{B}}\times\boldsymbol{\breve{g}}_\psi\boldsymbol{\cdot}\boldsymbol{\nabla} \breve{B}\; \frac{\iota - N/M}{{\left| \boldsymbol{\nabla}_{\varGamma} \alpha \right|}^2}\right) \right] , \end{gather}

(4.20b)\begin{gather}0=\int_{\mathcal{S}} \,\mathrm{d} S \left\{ q_\alpha \boldsymbol{\breve{B}} \boldsymbol{\cdot} \boldsymbol{\nabla}\phi + v_\mathrm{QS}\; \boldsymbol{\breve{B}}\times\boldsymbol{\breve{g}}_\psi\boldsymbol{\cdot}\boldsymbol{\nabla} \breve{B} \left[ \frac{\boldsymbol{\nabla}_{\varGamma}\alpha\boldsymbol{\cdot}\boldsymbol{\nabla}_{\varGamma} \phi}{{\left| \boldsymbol{\nabla}_{\varGamma} \alpha \right|}^2} (\iota-N/M) + 1 \right] \right\}. \end{gather}

Similarly to the rotational transform objective function case, $q_\alpha$ satisfies a magnetic differential equation (4.20a) on $\mathcal {S}$, with integral condition (4.20b). By (A 3), the surface integral of (4.20a) is consistent with the magnetic field's normal component vanishing on the boundary (4.2b).

Furthermore, requiring the Lagrangian to be stationary with respect to variations in $\omega$, we obtain the adjoint equations for $q_\omega$ from the third and fourth lines of (B 26):

(4.21a)\begin{equation} \Delta q_\omega = 0, \quad\mathrm{in}\ \mathcal{V}, \end{equation}

(4.21b)\begin{align} & \boldsymbol{\nabla} q_\omega \boldsymbol{\cdot} {\boldsymbol{\hat{{n}}}} ={-}\boldsymbol{\nabla}_{\varGamma} \boldsymbol{\cdot} \left\{ q_\alpha \boldsymbol{\nabla}_{\varGamma} \alpha + v_\mathrm{QS} \;\boldsymbol{\nabla}_{\varGamma} \breve{B} - \frac{\boldsymbol{\breve{B}}}{\breve{B}} \boldsymbol{\nabla}_{\varGamma} \boldsymbol{\cdot} (v_\mathrm{QS} \; \boldsymbol{\breve{B}})\right.\nonumber\\ & \quad \left.+ \left(\iota-N/M\right) \left[ v_\mathrm{QS} \; \boldsymbol{\breve{g}}_\psi\times\boldsymbol{\nabla}_{\varGamma}\breve{B} - \boldsymbol{\breve{B}}\; \boldsymbol{\nabla}_{\varGamma}\boldsymbol{\cdot}\left(\frac{1}{\breve{B}}v_\mathrm{QS} \boldsymbol{\breve{B}}\times\boldsymbol{\breve{g}}_\psi\right) \right] \vphantom{\left\{ q_\alpha \boldsymbol{\nabla}_{\varGamma} \alpha + v_\mathrm{QS} \;\boldsymbol{\nabla}_{\varGamma} \breve{B} - \frac{\boldsymbol{\breve{B}}}{\breve{B}} \boldsymbol{\nabla}_{\varGamma} \boldsymbol{\cdot} (v_\mathrm{QS} \; \boldsymbol{\breve{B}})\right.}\right\}, \quad \text{on } \mathcal{S}. \end{align}

Again, $q_\omega$ satisfies the Laplace equation in $\mathcal {V}$ (4.21a), with a normal boundary condition on $\mathcal {S}$ that is the tangential divergence of a vector tangential to the surface (4.21b). The boundary condition (4.21b) is consistent with the Laplace equation, as $\int _\mathcal {V} \,\mathrm {d} V\; \Delta q_\omega = \int _\mathcal {S} \,\mathrm {d} S\; \boldsymbol {\nabla } q_\omega \boldsymbol {\cdot } {\boldsymbol {\hat {{n}}}} = 0$, by (A 3).

Finally, we obtain the shape gradient from the last three lines of (B 26):

(4.22)\begin{align} \mathcal{G}_\mathrm{QS} & ={-} ({\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}\breve{B}) \boldsymbol{\nabla}_{\varGamma}\boldsymbol{\cdot}\left( v_\mathrm{QS} \boldsymbol{\breve{B}} \right) - v_\mathrm{QS}\;\left(\boldsymbol{\breve{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{\hat{{n}}}} - {\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\breve{B}} \right)\boldsymbol{\cdot}\boldsymbol{\nabla}_{\varGamma}\breve{B} \nonumber\\ & \quad + \left(\iota-N/M\right)\; {\left| \boldsymbol{\breve{g}}_\psi \right|} \; \boldsymbol{\breve{B}}\times\boldsymbol{\nabla}\breve{B}\boldsymbol{\cdot} \left[ {\left| \boldsymbol{\nabla}_{\varGamma}\alpha \right|} \boldsymbol{\nabla}_{\varGamma}\left( \frac{v_\mathrm{QS}}{{\left| \boldsymbol{\nabla}_{\varGamma}\alpha \right|}} \right)\right.\nonumber\\ & \quad + \left.{\boldsymbol{\hat{{n}}}}\; v_\mathrm{QS}\; \left( \frac{\boldsymbol{\nabla}_{\varGamma}\alpha\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\varGamma}\alpha}{{\left| \boldsymbol{\nabla}_{\varGamma} \alpha \right|}^2} -h \right) \right] \nonumber\\ & \quad + \boldsymbol{\breve{B}}\boldsymbol{\cdot}\boldsymbol{\nabla} q_\omega + q_\alpha \left( {\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\breve{B}} - \boldsymbol{\breve{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{\hat{{n}}}} \right)\boldsymbol{\cdot}\boldsymbol{\nabla}_{\varGamma}\alpha + \frac{h}{2} v_\mathrm{QS}^2, \end{align}

with $\delta \mathcal {L}_\mathrm {QS}[\delta \boldsymbol {x}] = \int _{\mathcal {S}} \,\mathrm {d} S \;( {\delta \boldsymbol {{x}}\boldsymbol {\cdot } {\boldsymbol {\hat {{n}}}}})\;\mathcal {G_\mathrm {QS}}$ and $h$ the summed curvature (A 4).

The shape gradient (4.22) for targeted quasi-helical symmetry with helicity $N/M = 5$ is shown in figure 3(a) for the example rotating ellipse case introduced in figure 1. The shape gradient obtained through adjoint methods is verified against a finite-difference evaluation in figure 3(b). The error is visibly small, as is attested by the small relative error of the shape gradient shown in figure 3(c). Convergence of the relative error for a parameter derivative in a random direction in $\varOmega$, evaluated with the adjoint method and with a forward finite-difference scheme, is shown in figure 3(d). Akin to the rotational transform figure of merit convergence study in figure 2(d), the error decreases linearly with $\epsilon _\mathrm {FD}$ until it plateaus due to finite Fourier or radial resolution. While the lowest resolution of $(N_\mathrm {tor}, M_\mathrm {pol}) = (8,8)$ seemed reasonable for the rotational transform figure of merit in figure 2(d), a higher resolution is clearly required for the quasi-symmetry figure of merit. This could be due to the fact that higher derivatives of the magnetic field are involved in the shape gradient for quasi-symmetry (4.22) than in that for rotational transform (4.16), through derivatives of $v_\mathrm {QS}$. The resulting fine-scale structure of $\mathcal {G}$ is harder to resolve with a truncated Fourier series. However, the relative errors in figures 2(d) and 3(d) are similarly small at the highest Fourier resolutions employed.

Figure 3. Shape gradient for the quasi-symmetry objective function with helicity $N/M=5$, evaluated through (a) adjoint methods and (b) a forward finite-difference scheme with step size $\epsilon _\mathrm {FD}=10^{-9}$, for the example rotating ellipse case introduced in figure 1, with Fourier resolution $(N_\mathrm {tor}, M_\mathrm {pol}) = (16,16)$. The relative error, defined as the absolute error normalised by the maximal absolute value of the adjoint shape gradient, is shown in (c). The convergence of the relative error in the parameter derivative (2.4) for a random direction in $\varOmega$ is shown in (d) as a function of the step size $\epsilon _\mathrm {FD}$ and Fourier resolution $(N_\mathrm {tor}, M_\mathrm {pol})$. The black dashed line indicates the linear scaling in $\epsilon _\mathrm {FD}$ expected from the employed forward finite-difference scheme.