2.1 Sources and constraints

To avoid unphysical constraints on $f_{1s}$ implied by the moment equations of (2.1) in the presence of a non-zero $E_{r}$ (Landreman et al. Reference Landreman, Smith, Mollén and Helander2014), particle and heat sources are added to the DKE (A 1),

(2.7) $$\begin{eqnarray}\displaystyle \mathbb{L}_{0s}\,f_{1s}-C_{s}(f_{1s})-f_{Ms}(x_{s}^{2}-{\textstyle \frac{5}{2}})S_{1s}^{f}(\unicode[STIX]{x1D713})-f_{Ms}(x_{s}^{2}-{\textstyle \frac{3}{2}})S_{2s}^{f}(\unicode[STIX]{x1D713})=\mathbb{S}_{0s}, & & \displaystyle\end{eqnarray}$$

where $S_{1s}^{f}(\unicode[STIX]{x1D713})$ and $S_{2s}^{f}(\unicode[STIX]{x1D713})$ are unknowns such that $S_{1s}^{f}$ provides a particle source and $S_{2s}^{f}$ provides a heat source. The collisionless trajectory operator in SFINCS coordinates is

(2.8) $$\begin{eqnarray}\displaystyle \mathbb{L}_{0s}=\dot{\boldsymbol{r}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}+{\dot{x}}_{s}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{s}}+\dot{\unicode[STIX]{x1D709}}_{s}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}_{s}}, & & \displaystyle\end{eqnarray}$$

and the inhomogeneous drive term is $\mathbb{S}_{0s}=-(\boldsymbol{v}_{\text{m}s}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713})\unicode[STIX]{x2202}f_{Ms}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}$ . The source functions are determined via the requirement that $\langle \int \text{d}^{3}v\,f_{1s}\rangle _{\unicode[STIX]{x1D713}}=0$ and $\langle \int \text{d}^{3}v\,x_{s}^{2}f_{1s}\rangle _{\unicode[STIX]{x1D713}}=0$ (i.e. $f_{1s}$ does not provide net density or pressure). So, the following system of equations is solved,

(2.9) $$\begin{eqnarray}\displaystyle \underbrace{\left[\begin{array}{@{}ccc@{}}\mathbb{L}_{0s}-C_{s} & -f_{Ms}(x_{s}^{2}-{\textstyle \frac{5}{2}}) & -f_{Ms}(x_{s}^{2}-{\textstyle \frac{3}{2}})\\ \mathbb{ L}_{1s} & 0 & 0\\ \mathbb{ L}_{2s} & 0 & 0\end{array}\right]}_{\mathbb{L}_{s}}\underbrace{\left[\begin{array}{@{}c@{}}f_{1s}\\ S_{1s}^{\,f}\\ S_{2s}^{\,f}\end{array}\right]}_{F_{s}}=\underbrace{\left[\begin{array}{@{}c@{}}\mathbb{S}_{0s}\\ 0\\ 0\end{array}\right]}_{\mathbb{S}_{s}}. & & \displaystyle\end{eqnarray}$$

The velocity-space averaging operations are denoted $\mathbb{L}_{1s}\,f_{1s}=\langle \int \text{d}^{3}v\,f_{1s}\rangle _{\unicode[STIX]{x1D713}}$ and $\mathbb{L}_{2s}\,f_{1s}=\langle \int \text{d}^{3}v\,f_{1s}x_{s}^{2}\rangle _{\unicode[STIX]{x1D713}}$ . The full multi-species system can be written as,

(2.10) $$\begin{eqnarray}\displaystyle \left[\begin{array}{@{}c@{}}\mathbb{L}_{1}\\ \vdots \\ \mathbb{L}_{N_{\text{species}}}\end{array}\right]\left[\begin{array}{@{}c@{}}F_{1}\\ \vdots \\ F_{N_{\text{species}}}\end{array}\right]=\left[\begin{array}{@{}c@{}}\mathbb{S}_{1}\\ \vdots \\ \mathbb{S}_{N_{\text{species}}}\end{array}\right]. & & \displaystyle\end{eqnarray}$$

Here the linear systems corresponding to each species as in (2.9) are coupled through the collision operator. We use the following notation to refer to the above system,

(2.11) $$\begin{eqnarray}\displaystyle \mathbb{L}F=\mathbb{S}. & & \displaystyle\end{eqnarray}$$

3.1 Continuous approach

Let $F=\{F_{s}\}_{s=1}^{N_{\text{species}}}$ be the set of unknowns computed with SFINCS before discretization, denoted by the column vector in (2.10) with $F_{s}$ given by (2.9). That is, $F$ consists of a set of $N_{\text{species}}$ distribution functions over $(\unicode[STIX]{x1D703},\unicode[STIX]{x1D701},x_{s},\unicode[STIX]{x1D709}_{s})$ and their associated source functions. We define an inner product between two such quantities in the following way,

(3.1) $$\begin{eqnarray}\displaystyle \langle F,G\rangle =\mathop{\sum }_{s}\left\langle \int \text{d}^{3}v\,\frac{f_{1s}g_{1s}}{f_{Ms}}\right\rangle _{\unicode[STIX]{x1D713}}+S_{1s}^{f}S_{1s}^{g}+S_{2s}^{f}S_{2s}^{g}. & & \displaystyle\end{eqnarray}$$

Here the superscript on $S_{1s}$ and $S_{2s}$ denotes the distribution function with which the source functions are associated and the sum is over species. The space of continuous functions, $F$ , of this form such that $\langle F,F\rangle$ is bounded will be denoted by ${\mathcal{H}}$ . It can be seen that (3.1) is indeed an inner product, as it satisfies conjugate symmetry ( $\langle G,F\rangle =\langle F,G\rangle$ $\forall F,G\in {\mathcal{H}}$ ), linearity ( $\langle F+G,H\rangle =\langle F,H\rangle +\langle G,H\rangle$ $\forall F,G,H\in {\mathcal{H}}$ and $\langle F,aG\rangle =a\langle F,G\rangle$ $\forall F,G\in {\mathcal{H}}$ , $a\in \mathbb{R}$ ), and positive definiteness ( $\langle F,F\rangle \geqslant 0$ and $\langle F,F\rangle =0$ only if $F=0$ $\forall F\in {\mathcal{H}}$ ) (Rudin Reference Rudin2006). This implies that if ${\mathcal{H}}$ is finite-dimensional, then for any linear operator $L$ there exists a unique adjoint operator $L^{\dagger }$ such that $\langle LF,G\rangle =\langle F,L^{\dagger }G\rangle$ for all $F,G\in {\mathcal{H}}$ . While here ${\mathcal{H}}$ is not finite-dimensional, we will show that such an adjoint operator exists for this inner product.

Note that the norm associated with this inner product $\Vert F\Vert =\sqrt{\langle F,F\rangle }$ is similar to the free energy norm,

(3.2) $$\begin{eqnarray}\displaystyle W=\mathop{\sum }_{s}\left\langle \int \text{d}^{3}v\,\frac{T_{s}\,f_{1s}^{2}}{2f_{Ms}}\right\rangle _{\unicode[STIX]{x1D713}}, & & \displaystyle\end{eqnarray}$$

which obeys a conservation equation in gyrokinetic theory (Krommes & Hu Reference Krommes and Hu1994; Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013; Landreman, Plunk & Dorland Reference Landreman, Plunk and Dorland2015). The choice of inner product (3.1) is advantageous, as the linearized Fokker–Planck collision operator becomes self-adjoint for species linearized about Maxwellians with the same temperature. In what follows, we assume that all included species are of the same temperature. This assumption could be lifted, with a modification to the collision operator that appears in the adjoint equation (see appendix B). This assumption is not necessary when using the discrete approach (see § 3.2).

Consider a moment of the distribution function ${\mathcal{R}}\in \{V_{\Vert ,s},\unicode[STIX]{x1D6E4}_{s},Q_{s},J_{b},J_{r},Q_{\text{tot}}\}$ , which can be written as an inner product with a vector $\widetilde{{\mathcal{R}}}\in {\mathcal{H}}$ ,

(3.3) $$\begin{eqnarray}\displaystyle {\mathcal{R}}=\langle F,\widetilde{{\mathcal{R}}}\rangle , & & \displaystyle\end{eqnarray}$$

according to (3.1). For example,

(3.4) $$\begin{eqnarray}\displaystyle \widetilde{J_{r}}=\left[\begin{array}{@{}c@{}}q_{s}\boldsymbol{v}_{\text{m}s}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}f_{Ms}\\ 0\\ 0\end{array}\right]_{s=1}^{N_{\text{species}}}, & & \displaystyle\end{eqnarray}$$

where the column structure corresponds with that in (2.9) and (2.10).

We are interested in computing the derivative of ${\mathcal{R}}$ with respect to a set of parameters, $\unicode[STIX]{x1D6FA}=\{\unicode[STIX]{x1D6FA}_{i}\}_{i=1}^{N_{\unicode[STIX]{x1D6FA}}}$ . This derivative can be computed with the chain rule,

(3.5) $$\begin{eqnarray}\displaystyle \left(\frac{\unicode[STIX]{x2202}{\mathcal{R}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}\right)_{\mathbb{L}F=\mathbb{S}}=\left(\frac{\unicode[STIX]{x2202}{\mathcal{R}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}\right)_{F}+\left\langle \widetilde{{\mathcal{R}}},\left(\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}\right)_{\mathbb{L}F=\mathbb{S}}\right\rangle . & & \displaystyle\end{eqnarray}$$

The subscripts in (3.5) denote the quantity that is held fixed while the derivative is computed. The first term on the right-hand side accounts for the explicit dependence on $\unicode[STIX]{x1D6FA}_{i}$ while the second accounts for the implicit dependence on $\unicode[STIX]{x1D6FA}_{i}$ through $F$ . Here $(\unicode[STIX]{x2202}F/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i})_{\mathbb{L}F=\mathbb{S}}$ can be computed by considering perturbations to the linear system (2.11), noting that in general both $\mathbb{L}$ and $\mathbb{S}$ can depend on $\unicode[STIX]{x1D6FA}$ ,

(3.6) $$\begin{eqnarray}\displaystyle \mathbb{L}\left(\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}\right)_{\mathbb{L}F=\mathbb{S}}=\left(\frac{\unicode[STIX]{x2202}\mathbb{S}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}-\frac{\unicode[STIX]{x2202}\mathbb{L}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}F\right). & & \displaystyle\end{eqnarray}$$

Computing $(\unicode[STIX]{x2202}F/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA})_{\mathbb{L}F=\mathbb{S}}$ using (3.6) requires solving $N_{\unicode[STIX]{x1D6FA}}$ linear systems of the same dimension as the DKE (2.11). To avoid this additional computational cost, we instead solve an adjoint equation,

(3.7) $$\begin{eqnarray}\displaystyle \mathbb{L}^{\dagger }q^{{\mathcal{R}}}=\widetilde{{\mathcal{R}}}. & & \displaystyle\end{eqnarray}$$

In what follows, we show that the adjoint variable, $q^{{\mathcal{R}}}$ , can be used to compute $(\unicode[STIX]{x2202}{\mathcal{R}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA})_{\mathbb{L}F=\mathbb{S}}$ without solving (3.6) for each $\unicode[STIX]{x1D6FA}_{i}$ . Using (3.7) with (3.5),

(3.8) $$\begin{eqnarray}\displaystyle \left(\frac{\unicode[STIX]{x2202}{\mathcal{R}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}\right)_{\mathbb{L}F=\mathbb{S}}=\left(\frac{\unicode[STIX]{x2202}{\mathcal{R}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}\right)_{F}+\left\langle \mathbb{L}^{\dagger }q^{{\mathcal{R}}},\left(\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}\right)_{\mathbb{L}F=\mathbb{S}}\right\rangle , & & \displaystyle\end{eqnarray}$$

and applying the adjoint property, we obtain

(3.9) $$\begin{eqnarray}\displaystyle \left(\frac{\unicode[STIX]{x2202}{\mathcal{R}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}\right)_{\mathbb{L}F=\mathbb{S}}=\left(\frac{\unicode[STIX]{x2202}{\mathcal{R}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}\right)_{F}+\left\langle q^{{\mathcal{R}}},\mathbb{L}\left(\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}\right)_{\mathbb{L}F=\mathbb{S}}\right\rangle . & & \displaystyle\end{eqnarray}$$

Using (3.6),

(3.10) $$\begin{eqnarray}\displaystyle \left(\frac{\unicode[STIX]{x2202}{\mathcal{R}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}\right)_{\mathbb{L}F=\mathbb{S}}=\left(\frac{\unicode[STIX]{x2202}{\mathcal{R}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}\right)_{F}+\left\langle q^{{\mathcal{R}}},\left(\frac{\unicode[STIX]{x2202}\mathbb{S}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}-\frac{\unicode[STIX]{x2202}\mathbb{L}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}F\right)\right\rangle . & & \displaystyle\end{eqnarray}$$

So, (3.10) provides the same derivative information as (3.5). Thus, using (3.10), the derivative with respect to $\unicode[STIX]{x1D6FA}$ can be computed with the solution to two linear systems, (2.11) and (3.7). The partial derivatives on the right-hand side of (3.10) can be computed analytically by considering the explicit geometric dependence of ${\mathcal{R}}$ , $\mathbb{L}$ and $\mathbb{S}$ .

When $N_{\unicode[STIX]{x1D6FA}}$ is large, the cost of computing $\unicode[STIX]{x2202}{\mathcal{R}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ using (3.10) is dominated not by the linear solve but by constructing $\unicode[STIX]{x2202}\mathbb{S}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ and $\unicode[STIX]{x2202}\mathbb{L}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ and computing the inner product. Thus the cost still scales with $N_{\unicode[STIX]{x1D6FA}}$ . However, we obtain a significant savings in comparison with forward difference derivatives, as shown in § 4.

The adjoint operator for each species takes the following form,

(3.11) $$\begin{eqnarray}\displaystyle \mathbb{L}_{s}^{\dagger }=\left[\begin{array}{@{}ccc@{}}\mathbb{L}_{0s}^{\dagger }-C_{s} & f_{Ms} & f_{Ms}x^{2}\\ \mathbb{L}_{1s}^{\dagger } & 0 & 0\\ \mathbb{L}_{2s}^{\dagger } & 0 & 0\end{array}\right], & & \displaystyle\end{eqnarray}$$

where $\mathbb{L}_{1s}^{\dagger }=5/2\mathbb{L}_{1s}-\mathbb{L}_{2s}$ and $\mathbb{L}_{2s}^{\dagger }=3/2\mathbb{L}_{1s}-\mathbb{L}_{2s}$ . The same column structure is used as for the forward operator (2.10), $\mathbb{L}^{\dagger }=\{\mathbb{L}_{s}^{\dagger }\}_{i=1}^{N_{\text{species}}}$ . The quantity $\mathbb{L}_{0s}^{\dagger }$ satisfies $\langle \int \text{d}^{3}v\,g_{1s}\mathbb{L}_{0s}\,f_{1s}/f_{Ms}\rangle _{\unicode[STIX]{x1D713}}=\langle \int \text{d}^{3}v\,f_{1s}\mathbb{L}_{0s}^{\dagger }g_{1s}/f_{Ms}\rangle _{\unicode[STIX]{x1D713}}$ and depends on which trajectory model is applied. The expression (3.11) can be verified by noting that

(3.12) $$\begin{eqnarray}\displaystyle \langle \mathbb{L}F,G\rangle & = & \displaystyle \mathop{\sum }_{s}\left\langle \frac{f_{1s}((\mathbb{L}_{0s}^{\dagger }-C_{s})g_{1s}+f_{Ms}(S_{1s}^{g}+S_{2s}^{g}x_{s}^{2}))}{f_{Ms}}\right\rangle _{\unicode[STIX]{x1D713}}+S_{1s}^{f}\mathbb{L}_{1s}^{\dagger }g_{1s}+S_{2s}^{f}\mathbb{L}_{2s}^{\dagger }g_{1s}\nonumber\\ \displaystyle & = & \displaystyle \langle F,\mathbb{L}^{\dagger }G\rangle .\end{eqnarray}$$

For the DKES trajectories the adjoint operator is,

(3.13) $$\begin{eqnarray}\displaystyle \mathbb{L}_{0s}^{\dagger }=-\mathbb{L}_{0s}. & & \displaystyle\end{eqnarray}$$

This anti-self-adjoint property is used in obtaining the variational principle which provides bounds on neoclassical transport coefficients in the DKES code (van Rij & Hirshman Reference van Rij and Hirshman1989). For full trajectories it is,

(3.14) $$\begin{eqnarray}\displaystyle \mathbb{L}_{0s}^{\dagger }=-\mathbb{L}_{0s}+\frac{q_{s}}{T_{s}}\frac{\text{d}\unicode[STIX]{x1D6F7}}{\text{d}\unicode[STIX]{x1D713}}\boldsymbol{v}_{\text{m}s}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}. & & \displaystyle\end{eqnarray}$$

The anti-self-adjoint property does not hold for this trajectory model as the $\boldsymbol{E}\times \boldsymbol{B}$ drift (C 9) is no longer divergenceless. See appendix C for details on obtaining these adjoint operators.

3.2 Discrete approach

Next we consider the discrete adjoint approach. Let $\overrightarrow{\boldsymbol{F}}$ be the set of unknowns computed with SFINCS after discretization of $F$ . The linear DKE (2.11) upon discretization can then be written schematically as

(3.15) $$\begin{eqnarray}\displaystyle \overleftrightarrow{\boldsymbol{L}}\overrightarrow{\boldsymbol{F}}=\overrightarrow{\boldsymbol{S}}. & & \displaystyle\end{eqnarray}$$

In this case, we can define an inner product as the vector dot product,

(3.16) $$\begin{eqnarray}\displaystyle \langle \overrightarrow{\boldsymbol{F}},\overrightarrow{\boldsymbol{G}}\rangle =\overrightarrow{\boldsymbol{F}}\boldsymbol{\cdot }\overrightarrow{\boldsymbol{G}}. & & \displaystyle\end{eqnarray}$$

In real Euclidean space, the adjoint operator, $(\overleftrightarrow{\boldsymbol{L}})^{\dagger }$ , which satisfies

(3.17) $$\begin{eqnarray}\displaystyle \langle \overleftrightarrow{\boldsymbol{L}}\overrightarrow{\boldsymbol{F}},\overrightarrow{\boldsymbol{G}}\rangle =\langle \overrightarrow{\boldsymbol{F}},(\overleftrightarrow{\boldsymbol{L}})^{\dagger }\overrightarrow{\boldsymbol{G}}\rangle & & \displaystyle\end{eqnarray}$$

is simply the transpose of the matrix, $(\overleftrightarrow{\boldsymbol{L}})^{\text{T}}$ . Again, the moments of the distribution function, ${\mathcal{R}}$ can be expressed as an inner product with a vector $\overrightarrow{\boldsymbol{R}}$ ,

(3.18) $$\begin{eqnarray}\displaystyle {\mathcal{R}}=\langle \overrightarrow{\boldsymbol{F}},\overrightarrow{\boldsymbol{R}}\rangle . & & \displaystyle\end{eqnarray}$$

Using the discrete approach, the following adjoint equation must be solved

(3.19) $$\begin{eqnarray}\displaystyle (\overleftrightarrow{\boldsymbol{L}})^{\text{T}}\overrightarrow{\boldsymbol{q}}^{{\mathcal{R}}}=\overrightarrow{\boldsymbol{R}}. & & \displaystyle\end{eqnarray}$$

The adjoint variable, $\overrightarrow{\boldsymbol{q}}^{{\mathcal{R}}}$ , can again be used to compute $(\unicode[STIX]{x2202}{\mathcal{R}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i})_{\overleftrightarrow{\boldsymbol{L}}\overrightarrow{\boldsymbol{F}}=\overrightarrow{\boldsymbol{S}}}$ ,

(3.20) $$\begin{eqnarray}\displaystyle \left(\frac{\unicode[STIX]{x2202}{\mathcal{R}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}\right)_{\overleftrightarrow{\boldsymbol{L}}\overrightarrow{\boldsymbol{F}}=\overrightarrow{\boldsymbol{S}}}=\left(\frac{\unicode[STIX]{x2202}{\mathcal{R}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}\right)_{\overrightarrow{\boldsymbol{F}}}+\left\langle \overrightarrow{\boldsymbol{q}}^{{\mathcal{R}}},\left(\frac{\unicode[STIX]{x2202}\overrightarrow{\boldsymbol{S}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}-\frac{\unicode[STIX]{x2202}\overleftrightarrow{\boldsymbol{L}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}\overrightarrow{\boldsymbol{F}}\right)\right\rangle . & & \displaystyle\end{eqnarray}$$

As with the continuous approach, the partial derivatives on the right-hand side can be computed analytically. In this way, the derivative of ${\mathcal{R}}$ with respect to $\unicode[STIX]{x1D6FA}$ can be computed with only two linear solves, (3.15) and (3.19).

In the SFINCS implementation, the DKE is typically solved with the preconditioned GMRES algorithm. In the continuous approach, a preconditioner matrix for both the forward and adjoint operator must be lower–upper ( $LU$ )-factorized. Here the preconditioner matrix is the same as the full matrix but without cross-species or speed coupling. As the adjoint matrix is sufficiently different from the forward matrix, we do not obtain convergence when the same preconditioner is used for both problems. However, in the discrete approach, the $LU$ -factorization for the preconditioner of the forward matrix can be reused for the preconditioner of the adjoint matrix (If a matrix $\unicode[STIX]{x1D63C}$ has been factorized as $\unicode[STIX]{x1D63C}=LU$ then $\unicode[STIX]{x1D63C}^{\text{T}}=U^{\text{T}}L^{\text{T}}$ where $U^{\text{T}}$ is lower triangular and $L^{\text{T}}$ is upper triangular). This provides a significant reduction in memory and computational cost for the discrete approach.

Furthermore, the discrete adjoint approach provides the exact derivatives for the discretized problem. With this method the adjoint equation is obtained using the vector dot product and matrix transpose which can be computed without any numerical approximation. The error in the derivatives obtained by the adjoint method are therefore only limited by the tolerance to which the linear solve is performed with GMRES. On the other hand, the continuous adjoint approach relies on a continuous inner product which must ultimately be approximated numerically. Thus the continuous approach provides the exact derivatives only in the limit that the discrete approximation of the inner product exactly reproduces the continuous inner product. Therefore we expect the results of the discrete and adjoint approaches to agree within discretization error, as will be demonstrated in § 4.

The continuous approach can be advantageous in that an adjoint equation may be prescribed independently of the discretization scheme. Note that in the discrete approach, the adjoint operator is obtained from the matrix transpose of the discretized forward operator, which implies that the same spatial and velocity resolution parameters must be used for both the forward and adjoint solutions. In this work we will employ the same discretization parameters for both the adjoint and forward problems, but this restriction is not required for the continuous approach.

5.1 Local magnetic sensitivity analysis

Using the adjoint method, it is possible to compute derivatives of a moment of the distribution function with respect to the Fourier amplitudes of the field strength, $\{\unicode[STIX]{x2202}{\mathcal{R}}/\unicode[STIX]{x2202}B_{mn}^{c}\}$ . Rather than consider sensitivity in Fourier space, we would like to compute the sensitivity to local perturbations of the field strength. We now quantify the relationship between these two representations of sensitivity information.

Consider the Gâteaux functional derivative (Delfour & Zolâsio Reference Delfour and Zolâsio2011) of ${\mathcal{R}}$ with respect to $B$ ,

(5.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}{\mathcal{R}}(\unicode[STIX]{x1D6FF}B;B(\boldsymbol{r}))=\lim _{\unicode[STIX]{x1D716}\rightarrow 0}\frac{{\mathcal{R}}(B(\boldsymbol{r})+\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FF}B(\boldsymbol{r}))-{\mathcal{R}}(B(\boldsymbol{r}))}{\unicode[STIX]{x1D716}}. & & \displaystyle\end{eqnarray}$$

Here we consider a perturbation to the field strength at fixed $I(\unicode[STIX]{x1D713})$ , $G(\unicode[STIX]{x1D713})$ and $\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})$ . As $\unicode[STIX]{x1D6FF}{\mathcal{R}}(\unicode[STIX]{x1D6FF}B;B(\boldsymbol{r}))$ is a linear functional of $\unicode[STIX]{x1D6FF}B$ , by the Riesz representation theorem (Rudin Reference Rudin2006), $\unicode[STIX]{x1D6FF}{\mathcal{R}}$ can be expressed as an inner product with $\unicode[STIX]{x1D6FF}B$ and some element of the appropriate space. The function $\unicode[STIX]{x1D6FF}B$ is defined on a flux surface, $\unicode[STIX]{x1D713}$ ; thus it is sensible to express $\unicode[STIX]{x1D6FF}{\mathcal{R}}$ in the following way,

(5.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}{\mathcal{R}}(\unicode[STIX]{x1D6FF}B;B(\boldsymbol{r}))=\langle S_{{\mathcal{R}}}\unicode[STIX]{x1D6FF}B(\boldsymbol{r})\rangle _{\unicode[STIX]{x1D713}}. & & \displaystyle\end{eqnarray}$$

Here $\unicode[STIX]{x1D6FF}B(\boldsymbol{r})$ describes the local perturbation to the field strength, and $\unicode[STIX]{x1D6FF}{\mathcal{R}}$ quantifies the corresponding change to the moment ${\mathcal{R}}$ . The function $S_{{\mathcal{R}}}$ is analogous to the shape gradient, which quantifies the change in a figure of merit which results from a differential perturbation to a shape (Landreman & Paul Reference Landreman and Paul2018). The shape gradient will be discussed further in § 5.2.2.

Suppose that $B$ is stellarator symmetric and $N_{P}$ symmetric. If $E_{r}=0$ , then $S_{{\mathcal{R}}}$ must also possess stellarator and $N_{P}$ symmetry (see appendix D). However, when $E_{r}\neq 0$ , $S_{{\mathcal{R}}}$ is no longer guaranteed to have stellarator symmetry. Nonetheless, it may be desirable to ignore the stellarator-asymmetric part of $S_{{\mathcal{R}}}$ if an optimized stellarator-symmetric configuration is desired. For the remainder of this work, we will make this assumption, though the analysis could be extended to consider the effect of breaking of stellarator symmetry. The quantity $S_{{\mathcal{R}}}$ can be approximated by a truncated Fourier series under these assumptions,

(5.3) $$\begin{eqnarray}\displaystyle S_{{\mathcal{R}}}=\mathop{\sum }_{k}S_{m_{k}n_{k}}\cos (m_{k}\unicode[STIX]{x1D703}-n_{k}\unicode[STIX]{x1D701}), & & \displaystyle\end{eqnarray}$$

where $k$ sums over $m\leqslant m_{\max }$ and $|n|\leqslant n_{\max }$ such that $n$ is an integer multiple of $N_{P}$ . The quantity $\unicode[STIX]{x1D6FF}B(\boldsymbol{r})$ can be written in terms of perturbations to the Fourier coefficients,

(5.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}B(\boldsymbol{r})=\mathop{\sum }_{j}\unicode[STIX]{x1D6FF}B_{m_{j}n_{j}}^{c}\cos (m_{j}\unicode[STIX]{x1D703}-n_{j}\unicode[STIX]{x1D701}), & & \displaystyle\end{eqnarray}$$

where again the sum is only taken over $N_{P}$ symmetric modes. Now $\unicode[STIX]{x1D6FF}{\mathcal{R}}$ can be written in terms of perturbations to the Fourier coefficients,

(5.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}{\mathcal{R}}=\mathop{\sum }_{j}\frac{\unicode[STIX]{x2202}{\mathcal{R}}}{\unicode[STIX]{x2202}B_{m_{j}n_{j}}^{c}}\unicode[STIX]{x1D6FF}B_{m_{j}n_{j}}^{c}. & & \displaystyle\end{eqnarray}$$

In this way, (5.2) can be expressed as a linear system,

(5.6) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}{\mathcal{R}}}{\unicode[STIX]{x2202}B_{m_{j}n_{j}}^{c}}=\mathop{\sum }_{k}D_{jk}S_{m_{k}n_{k}}, & & \displaystyle\end{eqnarray}$$

where

(5.7) $$\begin{eqnarray}\displaystyle D_{jk}=V^{\prime }(\unicode[STIX]{x1D713})^{-1}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D703}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D701}\sqrt{g}\cos (m_{j}\unicode[STIX]{x1D703}-n_{j}\unicode[STIX]{x1D701})\cos (m_{k}\unicode[STIX]{x1D703}-n_{k}\unicode[STIX]{x1D701}). & & \displaystyle\end{eqnarray}$$

If the same number of modes is used to discretize $\unicode[STIX]{x1D6FF}{\mathcal{R}}$ and $S_{{\mathcal{R}}}$ , then the linear system is square.

Figure 3. (a) The local magnetic sensitivity function for the bootstrap current, $S_{J_{b}}$ , is shown for the W7-X standard configuration. Positive values indicate that increasing the field strength at a given location will increase $J_{b}$ through (5.2). (b) The local sensitivity function for the ion particle flux, $S_{\unicode[STIX]{x1D6E4}_{i}}$ .

In contrast to derivatives with respect to the Fourier modes of $B$ , the sensitivity function, $S_{{\mathcal{R}}}$ , is a spatially local quantity, quantifying the change in a figure of merit resulting from a local perturbation of the field strength. In this way, $S_{{\mathcal{R}}}$ can inform where perturbations to the magnetic field strength can be tolerated. The sensitivity function could be related directly to a local magnetic tolerance using the method described in section 9 of Landreman & Paul (Reference Landreman and Paul2018). In contrast with that work, here we are considering perturbations to the field strength on any flux surface rather than at the plasma boundary. However, $S_{{\mathcal{R}}}$ still provides insight into where trim coils should be placed or coil displacements can be tolerated without sacrificing desired neoclassical properties. The sensitivity function can also be used for gradient-based optimization in the space of the field strength on a flux surface as in § 5.2.1.

We compute $S_{J_{b}}$ for the W7-X standard configuration at $\unicode[STIX]{x1D70C}=0.70$ , shown in figure 3(a). We use a fixed-boundary equilibrium that preceded the coil design and does not include coil ripple, and the full equilibrium is used rather than the truncated Fourier series considered in § 4. The same resolution parameters are used as in § 4, and derivatives with respect to $B_{mn}^{c}$ are computed for $m_{\max }=n_{\max }=20$ . The largest modes for this configuration are the helical curvature $B_{11}^{c}$ , the toroidal curvature $B_{10}^{c}$ and the toroidal mirror $B_{01}^{c}$ . We find that $S_{J_{b}}$ is large and negative on the inboard side, indicating that increasing the magnitude of the toroidal curvature component of $B$ would lead to an increase in $J_{b}$ . This result is in agreement with previous analysis of the dependence of the bootstrap current on these three modes in the W7-X magnetic configuration space (Maassberg, Lotz & Nührenberg Reference Maassberg, Lotz and Nührenberg1993), which found that at low collisionality the bootstrap current coefficients depend strongly on the toroidal curvature. In figure 3(b) is the sensitivity function for the ion particle flux, $S_{\unicode[STIX]{x1D6E4}_{i}}$ , computed for the same configuration using $m_{\max }=20$ and $n_{\max }=25$ . We find that the particle flux is more sensitive to perturbations on the outboard side in localized regions, while on the inboard side the sensitivity is relatively small in magnitude.

5.2 Gradient-based optimization

5.2.1 Optimization of field strength

As a second demonstration of the adjoint neoclassical method, we consider optimizing in the space of the field strength on a surface, taking $\unicode[STIX]{x1D6FA}=\{B_{mn}^{c}\}$ . As Boozer coordinates are used, the covariant form (4.1) satisfies $(\unicode[STIX]{x1D735}\times \boldsymbol{B})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}=0$ and the contravariant form (4.2) satisfies $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}=0$ . As we will artificially modify the field strength while keeping other geometry parameters fixed, the resulting field will not necessarily satisfy both of these conditions with both the covariant and contravariant forms. While there is no guarantee that the resulting field strength will be consistent with a global equilibrium solution, it provides insight into how local changes to the field strength can impact neoclassical properties. As a second step, the outer boundary could be optimized to match the desired field strength on a single surface. In § 5.2.2, we discuss how the derivatives computed in this work could be coupled to optimization of an MHD equilibrium.

We perform optimization with a Broyden–Fletcher–Goldfarb–Shanno (BFGS) quasi-Newton method (Nocedal & Wright Reference Nocedal and Wright1999) using an objective function $\unicode[STIX]{x1D712}^{2}=J_{b}^{2}$ . A backtracking line search is used at each iteration to find a step size that satisfies a condition of sufficient decrease of $\unicode[STIX]{x1D712}^{2}$ . We use the same equilibrium as in § 5.1, retaining modes $m\leqslant 12$ and $|n|\leqslant 12$ , and compute derivatives with respect to these modes. Convergence to $\unicode[STIX]{x1D712}^{2}\leqslant 10^{-10}$ was obtained within 8 BFGS iterations (28 function evaluations), as shown in figure 4(a). The difference in field strength between the initial and optimized configuration, $B_{\text{opt}}-B_{\text{init}}$ , is shown in figure 4(b). As expected from the analysis in § 5.1, the field strength increased on the outboard side and decreased on the inboard side in comparison with $B_{\text{init}}$ .

Figure 4. (a) Convergence of $\unicode[STIX]{x1D712}^{2}=J_{b}^{2}$ for optimization over $\unicode[STIX]{x1D6FA}=\{B_{mn}^{c}\}$ with an adjoint-based BFGS method. (b) The change in field strength from the initial to optimized configuration.

5.2.2 Optimization of MHD equilibria

The local sensitivity function, $S_{{\mathcal{R}}}$ , along with $\unicode[STIX]{x2202}{\mathcal{R}}/\unicode[STIX]{x2202}I$ , $\unicode[STIX]{x2202}{\mathcal{R}}/\unicode[STIX]{x2202}G$ , and $\unicode[STIX]{x2202}{\mathcal{R}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D704}$ , can be used to determine how perturbations to the outer boundary of the plasma, $\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}$ , result in perturbations to ${\mathcal{R}}$ . This is quantified through the idea of the shape gradient, which is described below. The partial derivatives of ${\mathcal{R}}$ can be computed with the adjoint method outlined in § 3, and the shape gradient can be obtained with only one additional MHD equilibrium solution through the application of another adjoint method.

Consider a figure of merit which is integrated over a toroidal domain, $\unicode[STIX]{x1D6E4}$ ,

(5.8) $$\begin{eqnarray}\displaystyle f_{{\mathcal{R}}}(\unicode[STIX]{x1D6E4})=\int _{\unicode[STIX]{x1D6E4}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713}){\mathcal{R}}(\unicode[STIX]{x1D713}), & & \displaystyle\end{eqnarray}$$

where $w(\unicode[STIX]{x1D713})$ is a weighting function. That is, SFINCS is run on a set of $\unicode[STIX]{x1D713}$ surfaces within $\unicode[STIX]{x1D6E4}$ and the volume integral is computed numerically. Here we consider $\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}$ to be the plasma boundary used for a fixed-boundary MHD equilibrium calculation. The perturbation to $f_{{\mathcal{R}}}$ resulting from normal perturbation to $\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}$ can be written in the following form,

(5.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}f_{{\mathcal{R}}}(\unicode[STIX]{x1D6E4};\unicode[STIX]{x1D6FF}\boldsymbol{r})=\int _{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}}\text{d}^{3}x\,(\unicode[STIX]{x1D6FF}\boldsymbol{r}\boldsymbol{\cdot }\boldsymbol{n}){\mathcal{G}}, & & \displaystyle\end{eqnarray}$$

under certain assumptions of smoothness (Delfour & Zolésio Reference Delfour and Zolésio2011). This can be thought of as another instance of the Riesz representation theorem, as $\unicode[STIX]{x1D6FF}f_{{\mathcal{R}}}$ is a linear functional of $\unicode[STIX]{x1D6FF}\boldsymbol{r}$ . Here $\boldsymbol{n}$ is the outward unit normal on $\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}$ and $\unicode[STIX]{x1D6FF}\boldsymbol{r}$ is a vector field describing the perturbation to the surface. Intuitively, only normal perturbations to $\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}$ result in a change to $f_{{\mathcal{R}}}$ . The shape gradient is ${\mathcal{G}}$ , which quantifies the contribution of a local normal perturbation of the boundary to the change in $f_{{\mathcal{R}}}$ . The shape gradient can be used for fixed-boundary optimization of equilibria or for analysis of sensitivity to perturbations of magnetic surfaces. It can be computed using a second adjoint method, where a perturbed MHD force balance equation is solved with the addition of a bulk force which depends on derivatives computed from the neoclassical adjoint method (Antonsen et al. Reference Antonsen, Paul and Landreman2019). While the continuous neoclassical adjoint method described in this work arises from the self-adjointness of the linearized Fokker–Planck operator, the adjoint method for MHD equilibria arises from the self-adjointness of the MHD force operator. In practice these two adjoint methods could be coupled by first computing an MHD equilibrium solution, computing neoclassical transport and its geometric derivatives from this equilibrium with the neoclassical adjoint method, and passing these derivatives back to the equilibrium code to compute the shape gradient with the perturbed MHD adjoint method. In this way derivatives of neoclassical quantities with respect to the shape of the outer boundary are computed with only two equilibrium solutions and two DKE solutions. This calculation will be reported in a future publication.

Rather than solve an additional adjoint equation, the outer boundary could be optimized by numerically computing derivatives of $\{B_{mn}^{c}(\unicode[STIX]{x1D713}),G(\unicode[STIX]{x1D713}),I(\unicode[STIX]{x1D713})\}$ with respect to the double Fourier series describing the outer boundary shape in cylindrical coordinates, $\{R_{mn}^{c},Z_{mn}^{s}\}$ , using a finite difference method. This could be done using the STELLOPT code (Spong et al. Reference Spong, Hirshman, Whitson, Batchelor, Carreras, Lynch and Rome1998; Reiman et al. Reference Reiman, Fu, Hirshman, Ku, Monticello, Mynick, Redi, Spong, Zarnstorff and Blackwell1999) with BOOZ_XFORM (Sanchez et al. Reference Sanchez, Hirshman, Ware, Berry and Spong2000) to perform the coordinate transformation. For example, if the rotational transform is held fixed in the variational moments equilibrium code (VMEC) equilibrium calculation (Hirshman & Whitson Reference Hirshman and Whitson1983), the derivative of a moment, ${\mathcal{R}}$ , with respect to a boundary coefficient, $R_{mn}^{c}$ , can be computed as,

(5.10) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}R_{mn}^{c}}=\mathop{\sum }_{m^{\prime }n^{\prime }}\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}B_{m^{\prime }n^{\prime }}^{c}(\unicode[STIX]{x1D713})}\frac{\unicode[STIX]{x2202}B_{m^{\prime }n^{\prime }}^{c}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}R_{mn}^{c}}+\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})}\frac{\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}R_{mn}^{c}}+\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}I(\unicode[STIX]{x1D713})}\frac{\unicode[STIX]{x2202}I(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}R_{mn}^{c}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})/\unicode[STIX]{x2202}B_{mn}^{c}(\unicode[STIX]{x1D713})$ , $\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})/\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})$ and $\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})/\unicode[STIX]{x2202}I(\unicode[STIX]{x1D713})$ are computed with the neoclassical adjoint method and $\unicode[STIX]{x2202}B_{mn}^{c}(\unicode[STIX]{x1D713})/\unicode[STIX]{x2202}R_{mn}^{c}$ , $\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})/\unicode[STIX]{x2202}R_{mn}^{c}$ and $\unicode[STIX]{x2202}I(\unicode[STIX]{x1D713})/\unicode[STIX]{x2202}R_{mn}^{c}$ are computed with finite difference derivatives using STELLOPT. Similarly, derivatives of $\{B_{mn}^{c}(\unicode[STIX]{x1D713}),G(\unicode[STIX]{x1D713}),I(\unicode[STIX]{x1D713})\}$ could be computed with respect to coil parameters using a free-boundary equilibrium solution, allowing for direct optimization of neoclassical quantities with respect to coil shapes. The neoclassical calculation with SFINCS is typically significantly more expensive than the equilibrium calculation (for the geometry discussed in § 5.1 fixed-boundary VMEC took 54 s while SFINCS took 157 s on 4 processors of the NERSC Edison computer). As such, combining adjoint-based with finite difference derivatives can still result in a significant computational savings.

5.3 Ambipolarity

As stellarators are not intrinsically ambipolar, the radial electric field is not truly an independent parameter. The ambipolar $E_{r}$ must be obtained which satisfies the condition $J_{r}(E_{r})=0$ . The application of adjoint-based derivatives for computing the ambipolar solution is discussed in § 5.3.1. An adjoint method to compute derivatives with respect to geometric parameters at fixed ambipolarity is discussed in § 5.3.2.

5.3.1 Accelerating ambipolar solve

A nonlinear root-finding algorithm must be used to compute the ambipolar $E_{r}$ . This root finding can be accelerated with derivative information, such as with a Newton–Raphson method (Press et al. Reference Press, Teukolsky, Vetterling and Flannery2007). The derivative required, $\unicode[STIX]{x2202}J_{r}/\unicode[STIX]{x2202}E_{r}$ , can be computed with the discrete or continuous adjoint method as described in § 3 with the replacement $\unicode[STIX]{x1D6FA}_{i}\rightarrow E_{r}$ , considering ${\mathcal{R}}=J_{r}$ .

Figure 5. The ambipolar root is obtained with Brent, Newton–Raphson and Newton hybrid nonlinear root solvers. The derivatives obtained with the adjoint method provide better convergence properties for the Newton methods.

We implement three nonlinear root-finding methods: Brent’s method (Brent Reference Brent2013), the Newton–Raphson method and a hybrid between the bisection and Newton–Raphson methods (Press et al. Reference Press, Teukolsky, Vetterling and Flannery2007). Brent’s method guarantees at least linear convergence by combining quadratic interpolation with bisection and does not require derivatives. The Newton–Raphson method can provide quadratic convergence under certain assumptions but in general is not guaranteed to converge. If an iterate lies near a stationary point or a poor initial guess is given, the method can fail. For this reason we implement the hybrid method, which combines the possible quadratic convergence properties of Newton–Raphson with the guaranteed linear convergence of the bisection method. Both Brent’s method and the hybrid method require the root to be bracketed, and therefore may require additional function evaluations in order to obtain the bracket.

We compare these methods in figure 5, using the W7-X standard configuration considered in § 5.1 with the full trajectory model and the discrete adjoint approach, beginning with an initial guess of $E_{r}=-10~\text{kV}~\text{m}^{-1}$ with bounds at $E_{r}^{\min }=-100~\text{kV}~\text{m}^{-1}$ and $E_{r}^{\max }=100~\text{kV}~\text{m}^{-1}$ . The root is located at $E_{r}=-3.56~\text{kV}~\text{m}^{-1}$ . For this example, the hybrid and Newton methods had nearly identical convergence properties, though the Newton method is less expensive as it does not require $J_{r}$ to be evaluated at the bounds of the interval. To obtain the same tolerance, the Newton method provided a 14 % savings in wall clock time over Brent’s method.

In the above discussion we have made the assumption that there is only one stable root of interest. Of course a given configuration may possess several roots, especially if the ions and electrons are in different collisionality regimes (Hastings, Houlberg & Shaing Reference Hastings, Houlberg and Shaing1985). Multiple roots can be obtained by performing several root solves with different initial values and brackets, which could be trivially parallelized. Thus the adjoint method could still provide an acceleration in this more general case.

5.3.2 Derivatives at ambipolarity

The adjoint method described in § 3 assumes that $E_{r}$ is held constant when computing derivatives with respect to $\unicode[STIX]{x1D6FA}$ . However, $E_{r}$ cannot truly be determined independently from geometric quantities, as the ambipolar solution should be recomputed as the geometry is altered. It is therefore desirable to compute derivatives at fixed ambipolarity (fixed $J_{r}=0$ ) rather than at fixed $E_{r}$ . This is performed by solving an additional adjoint equation,

(5.11) $$\begin{eqnarray}\displaystyle \mathbb{L}^{\dagger }q^{J_{r}}=\widetilde{J_{r}}, & & \displaystyle\end{eqnarray}$$

in the continuous approach or

(5.12) $$\begin{eqnarray}\displaystyle (\overleftrightarrow{\boldsymbol{L}})^{\text{T}}\overrightarrow{\boldsymbol{q}}^{J_{r}}=\overrightarrow{\boldsymbol{J}_{r}}, & & \displaystyle\end{eqnarray}$$

in the discrete approach. Details are described in appendix E.

It should be noted that by computing derivatives at ambipolarity we assume that a given moment ${\mathcal{R}}$ is a differentiable function of the geometry at fixed $J_{r}=0$ . That is, this method cannot be applied to cases in which a stable root disappears as the geometry varies. As this will occur at a stationary point of $J_{r}(E_{r})$ , this situation could be avoided within an optimization loop by computing derivatives at constant $E_{r}$ rather than constant $J_{r}$ if $|\unicode[STIX]{x2202}J_{r}/\unicode[STIX]{x2202}E_{r}|$ falls below a given threshold at ambipolarity.

Figure 6. (a) The cost of computing the gradient $\unicode[STIX]{x2202}{\mathcal{R}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ at ambipolarity scales with $N_{\unicode[STIX]{x1D6FA}}$ , the number of parameters in $\unicode[STIX]{x1D6FA}$ . (b) The fractional difference between $\unicode[STIX]{x2202}{\mathcal{R}}/\unicode[STIX]{x2202}B_{00}^{c}$ at constant ambipolarity obtained with the adjoint method and with finite difference derivatives.

Although an additional adjoint solve is required, this method of computing derivatives at ambipolarity is advantageous as several linear solves are typically required to obtain the ambipolar root. A comparison of the computational cost between the adjoint method and forward difference method for derivatives at ambipolarity is shown in figure 6(a). Here the full trajectory model is used, and the result for both the discrete and continuous adjoint methods are shown. For the finite difference derivative, the ambipolar solve is performed with Brent’s method at each step in $\unicode[STIX]{x1D6FA}$ . As in figure 2(b), we find that for large $N_{\unicode[STIX]{x1D6FA}}$ the cost of the continuous and discrete approaches are essentially the same, as the cost is no longer dominated by the linear solve. When computing the derivatives at ambipolarity, both adjoint methods decrease the cost by a factor of approximately $200$ for large $N_{\unicode[STIX]{x1D6FA}}$ .

Figure 7. The sensitivity function for the ion particle flux, $S_{\unicode[STIX]{x1D6E4}_{i}}$ , is computed at (a) constant $E_{r}$ and (b) constant $J_{r}$ . Similarly, $S_{J_{b}}$ is computed at (c) constant $E_{r}$ and (d) constant $J_{r}$ .

In figure 6(b) we show a benchmark between derivatives at ambipolarity, $(\unicode[STIX]{x2202}{\mathcal{R}}/\unicode[STIX]{x2202}B_{00}^{c})_{J_{r}}$ , computed with the discrete adjoint method and with forward difference derivatives. For the forward difference method, the Newton solver is used to obtain the ambipolar $E_{r}$ as $B_{00}^{c}$ is varied. As the forward difference step size $\unicode[STIX]{x0394}B_{00}^{c}$ decreases, the fractional difference again decreases proportional to $\unicode[STIX]{x0394}B_{00}^{c}$ until it reaches a minimum when $\unicode[STIX]{x0394}B_{00}^{c}/B_{00}^{c}$ is approximately $10^{-4}$ . In comparison with figure 1, we see that the minimum fractional difference is slightly larger at fixed ambipolarity than at fixed $E_{r}$ , as the tolerance parameters associated with the Newton solver introduce an additional source of error to the forward difference approach.

In figure 7(a,b) we compare the sensitivity function for the particle flux, $S_{\unicode[STIX]{x1D6E4}_{i}}$ , computed using derivatives at constant $E_{r}$ with that computed at constant $J_{r}$ . Here derivatives are computed using the discrete adjoint method with full trajectories, and the sensitivity function is constructed as described in § 5.1. The configuration and numerical parameters are the same as described in § 5.1. At constant $J_{r}$ the large region of increased sensitivity on the outboard side that appears at constant $E_{r}$ remains, though the overall magnitude of the sensitivity decreases. Thus it may be important to account for the effect of the ambipolar $E_{r}$ when optimizing for radial transport. In figure 7(c,d) we perform the same comparison for $S_{J_{b}}$ , finding the derivatives at fixed $E_{r}$ and at fixed $J_{r}$ to be virtually identical. This is to be expected, as numerical calculations of neoclassical transport coefficients for W7-X have found that the bootstrap coefficients are much less sensitive to $E_{r}$ than those for the radial transport (figures 18 and 26 in Beidler et al. (Reference Beidler, Allmaier, Isaev, Kasilov, Kernbichler, Leitold, Maaßberg, Mikkelsen, Murakami and Schmidt2011)). Furthermore, the bootstrap current in the $1/\unicode[STIX]{x1D708}$ regime is independent of $E_{r}$ , and the finite-collisionality correction is small for optimized stellarators, such as W7-X (Helander, Parra & Newton Reference Helander, Parra and Newton2017). Therefore, the ambipolarity corrections to the derivatives are less important for $J_{b}$ than for the radial transport.