A fast, accurate and easy to implement Kapur–Rokhlin quadrature scheme for singular integrals in axisymmetric geometries

Abstract

Many applications in magnetic confinement fusion require the efficient calculation of surface integrals with singular integrands. The singularity subtraction approaches typically used to handle such singularities are complicated to implement and low-order accurate. In contrast, we demonstrate that the Kapur–Rokhlin quadrature scheme is well-suited for the logarithmically singular integrals encountered for a toroidally axisymmetric confinement system, is easy to implement and is high-order accurate. As an illustration, we show how to apply this quadrature scheme for the efficient and accurate calculation of the normal component of the magnetic field due to the plasma current on the plasma boundary, via the virtual-casing principle.

1. Introduction

Integral formulations and integral equations are effective and popular tools for magnetostatic and magnetohydrodynamic problems in magnetic confinement fusion (Shafranov & Zakharov 1972; Zakharov 1973; Freidberg, Grossmann & Haas 1976; Hirshman, van Rij & Merkel 1986; Hirshman & Neilson 1986; Merkel 1986; Chance 1997; Lazerson, Sakakibara & Suzuki 2013; Ludwig et al. 2006; Ludwig, Rodrigues & Bizarro 2013; Drevlak et al. 2018; O'Neil & Cerfon 2018; Malhotra et al. 2019a; Pustovitov & Chukashev 2021). They have intuitive physical interpretations (Shafranov & Zakharov 1972; Zakharov 1973; Hirshman & Neilson 1986; Lazerson et al. 2013; Hanson 2015; Pustovitov & Chukashev 2021), provide geometric flexibility (Merkel 1986; Hirshman et al. 1986; Chance 1997; O'Neil & Cerfon 2018; Malhotra et al. 2019a) and often reduce the dimension of the unknown quantities one solves for, thus reducing the number of unknowns (Merkel 1986; Hirshman et al. 1986; Chance 1997; O'Neil & Cerfon 2018; Malhotra et al. 2019a). However, there typically is a price to pay for these advantages. Integral formulations often involve singular integrands, which are subtle to handle numerically (Freidberg et al. 1976; Merkel 1986; Atkinson 1997; Chance 1997; Ludwig et al. 2006, 2013; Klöckner et al. 2013; Kress 2014; Landreman & Boozer 2016; Ricketson et al. 2016; Malhotra et al. 2019a). The numerical difficulty of integrating these singular integrands depends on the nature of the singularity, the distribution of sources and the relative location of the evaluation points (often known as target points or observation points) with respect to the sources. In this article, we focus on the common situation in which we are trying to evaluate layer potentials at the source locations. This is, for example, the standard situation when applying Green's identity (Freidberg et al. 1976; Hirshman et al. 1986; Merkel 1986; Pustovitov 2008; Lee et al. 2015; Malhotra et al. 2019b).

In the fusion community, the numerical difficulty due to the singularity of the integrand is usually addressed via the method of singularity subtraction (Freidberg et al. 1976; Merkel 1986; Chance 1997; Ludwig et al. 2006, 2013). The method is robust, but leads to a quadrature scheme with low-order convergence. Furthermore, it is complicated to implement, and the chances of making mistakes in the derivation of the quadrature scheme or its numerical implementation are high. The purpose of our work is to demonstrate that for the singular layer potential integrals encountered in axisymmetric confinement devices, which can be reduced to line integrals of singular periodic functions, the Kapur–Rokhlin quadrature scheme (Kapur & Rokhlin 1997) is as simple to implement as the trapezoidal rule and is a scheme with high-order convergence, leading to low error for few quadrature points. For non-axisymmetric applications, we recently presented an efficient high-order quadrature scheme based on a different approach (Malhotra et al. 2019b), and alternative schemes may also provide good performance (Bruno & Garza 2020; Wu & Martinsson 2021). However, for axisymmetric cases, none of these schemes reduce to as simple and efficient a method as the Kapur–Rokhlin approach we present here.

As discussed previously, there are many situations in the study of axisymmetric magnetic confinement fusion devices for which the simplicity and accuracy of the Kapur–Roklin scheme could be demonstrated. For this article, we choose to focus on the evaluation of single-layer and double-layer potentials, which are the layer potentials appearing in Green's identity, and which we define precisely in § 2. Our first numerical test is a numerical verification of an identity for the double-layer potential. Our second numerical test focuses on the single-layer potential, which we evaluate for an application of the virtual-casing principle, to calculate the normal component of the magnetic field due to the plasma current at points on the plasma boundary (Shafranov & Zakharov 1972; Zakharov 1973; Hanson 2015).

The structure of this article is as follows. We introduce our mathematical notation, layer potential representations and Kapur–Rokhlin quadrature for singular integrals in § 2. In § 3, we give a brief summary of the virtual-casing principle, discuss the mathematical difficulties associated with the numerical evaluation of the virtual-casing integral and describe our method for addressing these difficulties in toroidally axisymmetric geometries. We prove in § 4 that the integrands we consider in this article are logarithmically singular, and can therefore be integrated to high accuracy with the Kapur–Rokhlin quadrature scheme, which we present in § 5. We demonstrate the accuracy and high order of convergence of the scheme for an application of the virtual-casing principle and for the evaluation of a double-layer potential in § 6, and summarise our work in § 7.

2. Mathematical background

2.1. Description of toroidal volumes and surfaces

Throughout our discussion of toroidal geometries, we make use of the standard, right-handed cylindrical coordinates $(r, \phi, z)$. At a point with toroidal angle $\phi$, we write the orthonormal unit vectors as $\boldsymbol {e}_r(\phi )$, $\boldsymbol {e}_\phi (\phi )$ and $\boldsymbol {e}_z$. With this notation, we emphasise the fact that the radial and azimuthal unit vectors depend on the toroidal angle.

In this article, we focus on axisymmetric geometries, which means that we only consider surfaces and volumes of revolution. We take the $z$-axis as the axis of revolution and define a simple closed curve $\gamma$ in the $(r,z)$ plane. By rotating this curve about the $z$-axis through the toroidal angle $\phi \in [0, 2{\rm \pi} ]$, we obtain a closed surface of revolution $\varGamma$. Its interior $\varOmega$ is the corresponding volume of revolution. We refer to $\gamma$ as the generating curve of $\varGamma$. It is parameterised by a single variable $t$, which we assume has period $L$. We denote the components of $\gamma$ in the $(r,z)$ plane by $(r(t), z(t))$, and we identify a point $\boldsymbol {y} \in \varGamma$ by its toroidal revolution angle $\phi$ and its generating curve parameter $t$. Correspondingly, we often write $\boldsymbol {y} = \boldsymbol {y}(\phi, t)$ to stress this parameterisation. Moreover, we assume that $\gamma$ is a $C^1$ curve which does not intersect the $z$-axis, in the sense that the derivatives $r'(t)$ and $z'(t)$ are continuous on $[0,L]$ and there exists $R_{\rm min}>0$ for which $r(t) \geqslant R_{\rm min}$ on $[0,L]$. Finally, we assume that $\gamma$ is oriented so that the vector $\boldsymbol {n}(\boldsymbol {y}(\phi, t)) = (\partial \boldsymbol {y} / \partial \phi ) \times (\partial \boldsymbol {y} / \partial t) / \mathcal {J}(t)$ is the unit outward normal to $\varOmega$ at $\boldsymbol {y}$. The quantity $\mathcal {J}(t)=\left \| (\partial \boldsymbol {y} / \partial \phi ) \times (\partial \boldsymbol {y} / \partial t) \right \|$ is the Jacobian of the parameterisation.

2.2. Single-layer and double-layer potentials for axisymmetric geometries

Layer potentials are fundamental tools in representing solutions to the partial differential equations that arise in magnetostatic and magnetohydrodynamic calculations for magnetic confinement fusion (Merkel 1986; Chance 1997; Ludwig et al. 2006, 2013; Landreman & Boozer 2016; Drevlak et al. 2018). Given a surface $\varGamma$ and a free-space Green's function $(\boldsymbol {x}, \boldsymbol {y}) \mapsto G(\boldsymbol {x}, \boldsymbol {y})$ for a partial differential equation, the single-layer operator $\mathcal {S}$ and the double-layer operator $\mathcal {D}$ are defined by (Guenther & Lee 1996)

(2.1a,b)\begin{align} [\mathcal{S} \sigma](\boldsymbol{x}) = \iint_\varGamma G(\boldsymbol{x}, \boldsymbol{y}) \sigma(\boldsymbol{y}) \,\mathrm{d} \varGamma(\boldsymbol{y}) \quad \text{and} \quad [\mathcal{D} \sigma](\boldsymbol{x}) = \iint_\varGamma \frac{\partial G(\boldsymbol{x}, \boldsymbol{y})}{\partial \boldsymbol{n}(\boldsymbol{y})} \sigma(\boldsymbol{y}) \,\mathrm{d} \varGamma(\boldsymbol{y}), \end{align}

respectively. In the double-layer representation, the quantity $\partial G(\boldsymbol {x}, \boldsymbol {y}) / \partial \boldsymbol {n}(\boldsymbol {y}) = \boldsymbol {n}(\boldsymbol {y}) \boldsymbol {\cdot } \boldsymbol \nabla _{\boldsymbol {y}} G(\boldsymbol {x}, \boldsymbol {y})$ is the derivative of the Green's function in the outward normal direction at $\boldsymbol {y}$. The function $\sigma$ is called the density function in this representation. Functions expressible as single-layer and double-layer potentials automatically satisfy the partial differential equation associated with the Green's function everywhere except the boundary $\varGamma$.

The case of Laplace's equation in three dimensions is particularly prevalent in magnetostatic and magnetohydrodynamic settings (Merkel 1986; Chance 1997; Landreman & Boozer 2016; Drevlak et al. 2018). Here, the Green's function is

(2.2)\begin{equation} G(\boldsymbol{x}, \boldsymbol{y}) = \frac 1{4 {\rm \pi}\left\| \boldsymbol{x} - \boldsymbol{y} \right\|}, \end{equation}

and functions expressible as $\mathcal {S} \sigma$ or $\mathcal {D} \sigma$ are harmonic on $\mathbb {R}^3 \setminus \varGamma$. Assuming the density $\sigma$ is also axisymmetric in the sense that $\partial \sigma / \partial \phi = 0$, one can analytically compute the part of the surface integral over the revolution angle $\phi \in [0, 2{\rm \pi} ]$. The resulting single-layer and double-layer integrals then are one-dimensional line integrals, and the resulting Green's function in the integrand can be expressed in terms of complete elliptic integrals (Ludwig et al. 2006, 2013; Jardin 2010). We provide an explicit expression in § 4.1, as we treat in detail the application of our method to the calculation of the virtual-casing principle. At this point, we just highlight the fact that the singularity in the Green's function when $\boldsymbol {x}=\boldsymbol {y}$ requires the use of specialised quadrature when the target $\boldsymbol {x}$ is located on the surface $\varGamma$, or regularisation methods (Freidberg et al. 1976; Merkel 1986; Chance 1997; Landreman & Boozer 2016; Drevlak et al. 2018; Malhotra et al. 2019b). The purpose of this article is to show that for applications in axisymmetric geometries, the Kapur–Rokhlin quadrature scheme is simpler to implement than the known regularisation methods used in the magnetic confinement community, and leads to high-order convergence.

2.3. Kapur–Rokhlin quadrature

The Kapur–Rokhlin quadrature rules (Kapur & Rokhlin 1997) are a collection of high-order schemes for computing

(2.3)\begin{equation} \int_0^b f(t) \,\mathrm{d} t \end{equation}

when $f$ has an integrable singularity at the origin of the form $\log t$ or $t^\lambda$ for $\lambda >-1$. We focus on the specific Kapur–Rokhlin scheme for a logarithmic singularity of the form $f(t) = p(t) \log t + q(t)$, where $p$ and $q$ need not have known formulae, but are assumed sufficiently smooth.

2.3.1. A Kapur–Rokhlin quadrature scheme for non-periodic functions

The Kapur–Rokhlin quadrature rules are corrections to the trapezoidal rule. For the standard trapezoidal rule with equal spacing $h = b/M$, the quadrature nodes for a non-singular integrand would be $t_i = i h$ for $i=0, \dots, M$. However, we omit the quadrature node $t_0=0$ because the integrand $f$ is singular there. This yields the punctured trapezoidal rule

(2.4)\begin{equation} \int_0^b f(t) \,\mathrm{d} t \approx h \left[ \sum_{i=1}^{M-1} f_i + \frac 12 f_M \right], \end{equation}

where we have written the shorthand $f_i = f(t_i)$. The punctured trapezoidal rule is low-order accurate when $f$ is singular. For example, when $f$ is logarithmically singular at the origin, the punctured trapezoidal rule error typically decays, according to Martinsson (2019), as $O(|h\log h|)$.

The Kapur–Rokhlin corrections place additional quadrature nodes outside the integration domain $[0,b]$. Specifically, the corrections depend on a convergence rate parameter $n$ and a smoothness parameter $m$. One chooses $m \geqslant 3$ an odd integer subject to the constraint that $p$ and $q$ are $m$ times continuously differentiable on a wider interval $[-nh,b+mh]$. Under these conditions, the Kapur–Rokhlin scheme sets constants $\gamma _j$ for $j=\pm 1, \dots, \pm n$ and $\beta _l$ for $l = 1, \dots, (m-1)/2$ and defines the quadrature rule

(2.5)\begin{equation} T_{m,n}^M(f) = h \left[\sum_{i=1}^{M-1} f_i + \frac 12 f_M \right] + h \sum_{l=1}^{(m-1)/2} \beta_l (f_{M-l} - f_{M+l}) + h \sum_{1 \leqslant |j| \leqslant n} \gamma_j f_j \end{equation}

for $M \geqslant n + (m-1)/2$. The error obeys the asymptotic rate

(2.6)\begin{equation} \left| T_{m,n}^M(f) - \int_0^b f(t) \,\mathrm{d} t \right| = O(h^n) \end{equation}

as $h \to 0$. The endpoint correction coefficients $\{\gamma _j\}$ and $\{\beta _l\}$ are derived by analysing the Euler–Maclaurin formula for quadrature error and solving a linear system for correction coefficients to obtain high-order accuracy. Figure 1 illustrates this quadrature scheme and compares it with the punctured trapezoidal rule.

Figure 1. Nodes and weights for the Kapur–Rokhlin and punctured trapezoidal quadrature rules for estimating $\int _0^1 f(t) \,\mathrm {d} t$ with $f(t) = \cos (4{\rm \pi} t) \log |t| + t$. We have used $M=10$, $n=2$, and $m=9$. Note that $\beta _4 \approx -3 \times 10^{-4}$, so the weight corrections corresponding to $t_{6}=0.6$ and $t_{14}=1.4$ are not visually discernible.

2.3.2. A simplified quadrature for periodic functions

This subsection follows an argument nearly verbatim from Hao et al. (2014). We explain how the Kapur–Rokhlin scheme we just presented simplifies when computing

(2.7)\begin{equation} \int_{{-}b}^b f(t)\, \mathrm{d} t \end{equation}

when $f$ is $2b$-periodic and logarithmically singular at the origin. We may express these assumed properties of $f$ through the form

(2.8)\begin{equation} f(t) = p(t) \log \left|\sin \frac{{\rm \pi} t}{2b} \right| + q(t), \end{equation}

for $2b$-periodic functions $p$ and $q$.

We sum two applications of the original Kapur–Rokhlin scheme: one for

(2.9)\begin{equation} I_1 = \int_0^b f(t) \,\mathrm{d} t \end{equation}

and another for

(2.10)\begin{equation} I_2 = \int_{{-}b}^0 f(t) \,\mathrm{d} t = \int_0^b f({-}t) \,\mathrm{d} t. \end{equation}

We assume that $p, q \in C^m[-b, b]$ and generate $2M-1$ equispaced quadrature nodes with spacing $h=b/M$, defined by $t_i = ih$ for $i=\pm 1, \dots, \pm (M-1), M$. The corrected scheme for $I_1$ is

(2.11)\begin{equation} I_1 = h \left[\sum_{i=1}^{M-1} f_i + \frac 12 f_M \right] + h \sum_{l=1}^{(m-1)/2} \beta_l (f_{M-l} - f_{M+l}) + h \sum_{1 \leqslant |j| \leqslant n} \gamma_j f_j + O(h^n), \end{equation}

and the scheme for $I_2$ is

(2.12)\begin{equation} I_2 = h \left[\sum_{i=1}^{M-1} f_{{-}i} + \frac 12 f_{{-}M} \right] + h \sum_{l=1}^{(m-1)/2} \beta_l (f_{{-}M+l} - f_{{-}M-l}) + h \sum_{1 \leqslant |j| \leqslant n} \gamma_j f_{{-}j} + O(h^n). \end{equation}

By periodicity, we may identify $f_i$ with $f_{i+2M}$ for all $i$. It follows that an $n$th-order quadrature for $I_1 + I_2$ is

(2.13)\begin{align} I_1 + I_2 & = h \left[ \sum_{1 \leqslant |i| \leqslant M-1} f_i + f_M \right] + h \sum_{1 \leqslant |j| \leqslant n} \gamma_j (f_j + f_{{-}j})+ O(h^n) \nonumber\\ & = \sum_{\substack{j={-}M+1\nonumber\\ j \ne 0}}^M w_j f_j + O(h^n) \end{align}

with

(2.14)\begin{equation} w_j = \begin{cases} h(1+\gamma_j + \gamma_{{-}j}), & 1 \leqslant |j| \leqslant n, \\ h, & \text{otherwise}. \end{cases} \end{equation}

We note that, by periodicity, the endpoint corrections corresponding to the constants $\beta _l$ exactly cancel. Moreover, we no longer need to place additional quadrature nodes beyond $[-b,b]$ because periodicity identifies the new quadrature nodes with existing nodes.

Given a table of $\gamma _j$ weights, this quadrature scheme is as easy to implement as the trapezoidal rule and yields high-order convergence, with a quadrature error that is $O(h^n)$. The necessary $\gamma _j$ weights can be found in Kapur & Rokhlin (1997, table 6) for $n=2,6,10$. Figure 2 can be compared with figure 1 to view the simplifications for periodic integrands.

Figure 2. Nodes and weights for the Kapur–Rokhlin and punctured trapezoidal quadrature rules for estimating $\int _{t_0-{\rm \pi} }^{t_0+{\rm \pi} } f(t) \,\mathrm {d} t$. The function $f(t)$ is the $2{\rm \pi}$-periodic integrand in (6.5) for one of our later numerical tests with a logarithmic singularity at $t_0=1$. We have used $M=10$ and $n=2$.

3. The virtual-casing principle for toroidally axisymmetric domains

3.1. Formulation of the virtual-casing principle

For axisymmetric confinement devices, the virtual-casing principle is most often used to compute the poloidal flux or the poloidal magnetic field due to the toroidal current flowing in the plasma (Shafranov & Zakharov 1972; Zakharov 1973; Hirshman & Neilson 1986; Zakharov & Pletzer 1999). A poloidal magnetic field ${\boldsymbol {B}}^{{\rm pol}}$ at any point $\boldsymbol {y}(\phi, t) \in \varGamma$ can be expressed in terms of its poloidal flux function $\psi (r,z)$ and the parameterisation $(\phi,t) \mapsto \boldsymbol {y}(\phi, t)$ by (Freidberg 2014)

(3.1)\begin{align} \boldsymbol{B}^{{\rm pol}}(\boldsymbol{y}(\phi,t)) & = \boldsymbol\nabla \psi(r(t),z(t)) \times \boldsymbol\nabla \phi \nonumber\\ & = \boldsymbol\nabla \psi(r(t),z(t)) \times \left(\frac {\boldsymbol{e}_\phi(\phi)}{r(t)} \right) \nonumber\\ & ={-}\frac 1{r(t)} \dfrac{\partial \psi}{\partial z}(r(t), z(t)) \boldsymbol{e}_r(\phi) + \frac 1{r(t)} \dfrac{\partial \psi}{\partial r}(r(t), z(t)) \boldsymbol{e}_z. \end{align}

Consider an axisymmetric plasma confined by external coils in equilibrium. The poloidal field $\boldsymbol {B}^{{\rm pol}}$ at any location is the sum of the poloidal field $\boldsymbol {B}_{{\rm ext}}^{{\rm pol}}$ due to the external coils and of the poloidal field $\boldsymbol {B}_V^{{\rm pol}}$ due to the plasma current. The field $\boldsymbol {B}_V^{{\rm pol}}$ is given for all $\boldsymbol {x} \in \mathbb {R}^3$ by the Biot–Savart law:

(3.2)\begin{equation} \boldsymbol{B}_V^{\rm pol}(\boldsymbol{x})=\frac{\mu_{0}}{4{\rm \pi}} \iiint_{\varOmega}J^{{\rm tor}}_{V}(\boldsymbol{y}) \boldsymbol{e}_\phi(\boldsymbol{y}) \times \frac{\boldsymbol{x} - \boldsymbol{y}}{\left\| \boldsymbol{x} - \boldsymbol{y} \right\|^3} \,\mathrm{d} \boldsymbol{y}, \end{equation}

where $\mu _{0}$ is the permeability of free space and $J^{{\rm tor}}_{V}$ is the toroidal current density in the plasma. Equation (3.2) is a volume integral, which is expensive to evaluate numerically. The virtual-casing principle gives a formula for $\boldsymbol {B}^{{\rm pol}}_V$ that depends only on the full field $\boldsymbol {B}^{{\rm pol}}$ at the plasma boundary, and only requires the evaluation of a surface integral (i.e. line integral for axisymmetric domains) (Shafranov & Zakharov 1972; Zakharov 1973; Hanson 2015). In experimental settings, such a representation is useful because only the total magnetic field $\boldsymbol {B}^{{\rm pol}}$ may be directly measurable (Hutchinson 2002). In theoretical settings, the total magnetic field may be computed by solving the Grad–Shafranov equation (Grad & Rubin 1958; Shafranov 1958) numerically with a fixed-boundary solver. Specifically, the virtual-casing principle states that if $\varGamma$ is the flux surface bounding the plasma, then $\boldsymbol {B}_V^{{\rm pol}}$ can be written in terms of a field generated by the toroidal surface current $\boldsymbol {J}^{{\rm tor}}_S$ such that $\mu _0 \boldsymbol {J}^{{\rm tor}}_S = -\boldsymbol {n} \times \boldsymbol {B}^{{\rm pol}}$, according to Hanson (2015):

(3.3)\begin{equation} \boldsymbol{B}^{{\rm pol}}_V(\boldsymbol{x}) = \frac 1{4{\rm \pi}} \iint_{\varGamma} \left[ \frac{(\boldsymbol{n}(\boldsymbol{y}) \times \boldsymbol{B}^{{\rm pol}}(\boldsymbol{y})) \times (\boldsymbol{x} - \boldsymbol{y})}{\left\| \boldsymbol{x} - \boldsymbol{y} \right\|^3} \right] \mathrm{d} \varGamma(\boldsymbol{y}) + \begin{cases} \boldsymbol{B}^{{\rm pol}}(\boldsymbol{x}) & \boldsymbol{x} \in \varOmega \\ \boldsymbol{B}^{{\rm pol}}(\boldsymbol{x}) / 2 & \boldsymbol{x} \in \varGamma \\ \boldsymbol{0} & \boldsymbol{x} \notin \bar{\varOmega}. \end{cases} \end{equation}

For certain applications, one is only interested in the normal component of the poloidal magnetic field (Merkel 1986; Hirshman et al. 1986; Merkel 1987; Landreman 2017; Zhu et al. 2018). The previous equation then leads to the more compact form

(3.4)\begin{equation} \boldsymbol{n}(\boldsymbol{x}) \boldsymbol{\cdot} \boldsymbol{B}^{{\rm pol}}_V(\boldsymbol{x}) = \frac 1{4{\rm \pi}} \boldsymbol{n}(\boldsymbol{x}) \boldsymbol{\cdot} \iint_{\varGamma} \left[ \frac{(\boldsymbol{n}(\boldsymbol{y}) \times \boldsymbol{B}^{{\rm pol}}(\boldsymbol{y})) \times (\boldsymbol{x} - \boldsymbol{y})}{\left\| \boldsymbol{x} - \boldsymbol{y} \right\|^3} \right] \mathrm{d} \varGamma(\boldsymbol{y}) \end{equation}

for $\boldsymbol {x} \in \varGamma$.

The reduction of the integral necessary to compute the field or its normal component from a volume integral to a surface integral is convenient from the point of view of the limited number of values that need to be specified as inputs, and also from the point of view of the computational cost of the integration (Lazerson et al. 2013). The surface integral in (3.3) is significantly faster to evaluate than the volume integral (3.2), although certain codes still choose to compute the latter (Hanson et al. 2009; Marx & Lütjens 2017).

For axisymmetric situations, one may further take advantage of the axisymmetry of $\boldsymbol {B}^{{\rm pol}}$ to integrate with respect to $\phi$ analytically, and reduce (3.3) and (3.4) to one-dimensional integrals, which are even less computationally expensive. However, one encounters a mathematical and computational difficulty if one does so, because the surface integrals in (3.3) and (3.4) are in fact improper integrals, which must be understood in the Cauchy principal value sense. This is what we discuss in the following section.

3.2. Numerical evaluation of the normal component of the virtual-casing magnetic field in axisymmetric systems

3.2.1. Circumventing integrals in the principal value sense

Most applications in magnetic confinement fusion rely on version (3.4) of the virtual-casing principle, in which one wants to compute the normal component of the poloidal magnetic field $\boldsymbol {B}_V^{{\rm pol}}$ at a boundary point $\boldsymbol {x}\in \varGamma$. One could, in principle, calculate this normal component by first computing all the components of $\boldsymbol {B}_V^{{\rm pol}}$ by the virtual-casing principle (3.3), and then computing $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {B}_V^{{\rm pol}}$ by a straightforward inner product. In other words, one first evaluates the double integral in (3.4), and then takes the dot product with the normal vector $\boldsymbol {n}$ to the surface $\varGamma$ at the point $\boldsymbol {x}$ of interest. This method has the advantage that it produces all components of the poloidal magnetic field $\boldsymbol {B}_V^{{\rm pol}}$ as intermediary results. However, it also has the disadvantage that one must use a careful principal value integration procedure to interpret the virtual-casing principle for $\boldsymbol {x} \in \varGamma$, as discussed in Theorem A.3 and its proof in the Appendix. The high-order singularity cancellation quadrature scheme we recently proposed for singular integrals on general non-axisymmetric surfaces (Malhotra et al. 2019b) automatically yields the appropriate principal value of the integral. However, it does so thanks to the intrinsic two-dimensional nature of the integral. It does not reduce to a simple and efficient one-dimensional quadrature scheme for the principal value of the integral, as is needed for axisymmetric applications. Similarly, we are not aware of a version of the Kapur–Rokhlin quadrature scheme designed to calculate the Cauchy principal value of the virtual-casing integral.

To address this difficulty, we consider an alternative method to compute the normal poloidal field $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {B}_V^{{\rm pol}}$, based on calculating the vector potential $\boldsymbol {A}_S$ produced by the surface current $\boldsymbol {J}^{{\rm tor}}_S$, and then obtaining $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {B}_V^{\rm pol}$ as the tangential derivative of the poloidal flux $\psi _S$ caused by the surface currents, which is easily expressed in terms of $\boldsymbol {A}_S$. This method is also used, in a slightly different form, by the stellarator optimisation code ROSE (Drevlak et al. 2018), but ROSE does not combine it with a high-order quadrature scheme. The vector potential is defined by

(3.5)\begin{equation} \boldsymbol{A}_S(\boldsymbol{x}) = \frac {\mu_0}{4 {\rm \pi}} \iint_\varGamma \frac{\boldsymbol{J}^{{\rm tor}}_S(\boldsymbol{y})}{\left\| \boldsymbol{x} - \boldsymbol{y} \right\|} \,\mathrm{d} \varGamma(\boldsymbol{y}) \end{equation}

and we recognise this expression as a vector of component-wise single-layer potentials for Laplace's equation in three dimensions, as introduced in § 2.2.

In contrast to the first approach, our method does not produce all components of $\boldsymbol {B}^{{\rm pol}}_V$. On the other hand, it only requires weakly singular integrals, because it only requires evaluations of single-layer potentials in (3.5), as we prove in § 4. In the next section, we describe this alternative method for computing the normal magnetic field in detail.

3.2.2. Normal component of the magnetic field as a tangential derivative

As in the derivation of the virtual-casing principle (Shafranov & Zakharov 1972; Zakharov 1973; Hanson 2015), one can interpret the plasma boundary $\varGamma$ of the axisymmetric plasma as a perfectly conducting shell which contributes to confining the plasma via the poloidal magnetic field $\boldsymbol {B}^{{\rm pol}}_S$ generated by the toroidal surface current density $\boldsymbol {J}^{{\rm tor}}_S$ flowing in $\varGamma$. At equilibrium, the poloidal surface current field $\boldsymbol {B}^{{\rm pol}}_S$ matches the poloidal field $\boldsymbol {B}_{{\rm ext}}^{{\rm pol}}$ from the external confining coils. The decomposition $\boldsymbol {B}^{{\rm pol}} = \boldsymbol {B}^{{\rm pol}}_V + \boldsymbol {B}^{{\rm pol}}_S$ then immediately yields $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {B}^{{\rm pol}}_V = -\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {B}^{{\rm pol}}_S$. Now, recall that we may represent $\boldsymbol {B}^{\rm pol}_S$ in axisymmetry via (3.1) to express its normal component at $(r, \phi, z)$ as

(3.6)\begin{equation} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{B}^{{\rm pol}}_S = \boldsymbol{n} \boldsymbol{\cdot} (\boldsymbol\nabla \psi_S \times \boldsymbol\nabla \phi) ={-}\frac{1}{r} (\boldsymbol{n} \times \boldsymbol{e}_\phi(\phi)) \boldsymbol{\cdot} \boldsymbol\nabla \psi_S \end{equation}

by the circular shift identity of the vector triple product. Let $\boldsymbol {x}$ correspond to the parameter pair $(\phi _0, t_0)$ with $\boldsymbol {x} = (R, \phi _0, Z) = (r(t_0), \phi _0, z(t_0))$. We define the tangent vector $\boldsymbol {t}$ at a point $\boldsymbol {x} \in \varGamma$ as

(3.7)\begin{align} \boldsymbol{t}(\boldsymbol{x}) & = \boldsymbol{n}(\boldsymbol{x}) \times \boldsymbol{e}_\phi(\phi_0) \nonumber\\ & = \frac{(z'(t_0) \boldsymbol{e}_r(\phi_0) - r'(t_0) \boldsymbol{e}_z ) \times \boldsymbol{e}_\phi(\phi_0)} {\sqrt{r'(t_0)^2 + z'(t_0)^2}} \nonumber\\ & = \frac{ r'(t_0) \boldsymbol{e}_r(\phi_0) + z'(t_0) \boldsymbol{e}_z } {\sqrt{r'(t_0)^2 + z'(t_0)^2}}. \end{align}

It follows that

(3.8)\begin{align} \boldsymbol{n}(\boldsymbol{x}) \boldsymbol{\cdot} \boldsymbol{B}^{\rm pol}_V(\boldsymbol{x}) & = \frac{1}{r(t_0)} \boldsymbol{t}(\boldsymbol{x}) \boldsymbol{\cdot} \boldsymbol\nabla \psi_S(r(t_0), z(t_0)) \nonumber\\ & = \frac{1}{r(t_0)\sqrt{r'(t_0)^2 + z'(t_0)^2}} \left( \frac{\partial{\psi_S}}{\partial{r}} r'(t_0) + \frac{\partial{\psi_S}}{\partial{z}} z'(t_0) \right)\nonumber\\ & = \frac{1}{\mathcal{J}(t_0)} \frac{\mathrm{d} \psi_S(r(t_0), z(t_0))}{\mathrm{d} t_0}, \end{align}

where $\mathcal {J}(t_0) = r(t_0)\sqrt {r'(t_0)^2 + z'(t_0)^2}$ is the Jacobian of the parameterisation introduced in § 2.1.

Finally, the poloidal flux function and the vector potential are related at $\boldsymbol {x}$ by the relation (Jardin 2010; Freidberg 2014)

(3.9)\begin{equation} \psi_S(R, Z) = R \boldsymbol{e}_\phi(\phi_0) \boldsymbol{\cdot} \boldsymbol{A}_S(\boldsymbol{x}). \end{equation}

Given only point evaluations of $\boldsymbol {A}_S$ at equispaced parameters $\{t_i\}$, one can compute $\mathrm {d} \psi _S / \mathrm {d} t$ with high-order accuracy at each $t_i$ by Fourier differentiation. This leads to our high-order accurate approach to the virtual-casing principle in axisymmetric geometries. One computes $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {B}^{{\rm pol}}_V$ on $\varGamma$ by evaluating a weakly singular integral with a high-order accurate quadrature rule, and then Fourier differentiating. We show in § 5 that Kapur–Rokhlin quadrature provides a simple way to obtain high-order accuracy for the weakly singular integral. Before we do so, we prove in the next section that it indeed is a weakly singular integral, whose properties satisfy the requirements to obtain high accuracy with the Kapur–Rokhlin quadrature rule.

4. Analytical results for singular integrals

4.1. Vector potential singularity

By the virtual-casing principle, we have $\mu _0 \boldsymbol {J}^{{\rm tor}}_S = - \boldsymbol {n} \times \boldsymbol {B}^{{\rm pol}}$, so we may equivalently write the vector potential as

(4.1)\begin{equation} \boldsymbol{A}_S(\boldsymbol{x}) ={-}\frac{1}{4{\rm \pi}} \iint_\varGamma \frac{\boldsymbol{n}(\boldsymbol{y}) \times \boldsymbol{B}^{{\rm pol}}(\boldsymbol{y})} {\left\| \boldsymbol{x} - \boldsymbol{y} \right\|} \mathrm{d} \varGamma(\boldsymbol{y}), \end{equation}

which satisfies $\boldsymbol {B}^{{\rm pol}}_S = \boldsymbol { \nabla } \times \boldsymbol {A}_S$. As mentioned earlier, we may view this expression as a vector of component-wise single-layer potentials for Laplace's equation in three dimensions. Physically, it is clear that the only non-zero component of the vector potential is in the $\boldsymbol {e}_\phi$ direction. Furthermore, Chance (1997, § V) shows that after integrating in the toroidal angle $\phi$, the one-dimensional integrands of both the single-layer potential $\mathcal {S}[\mu _0\boldsymbol {J}^{{\rm tor}}_S]$ and the double-layer potential $\mathcal {D}[\mu _0\boldsymbol {J}^{{\rm tor}}_S]$ have integrable, logarithmic singularities. Chance shows this in a modified coordinate system using the poloidal flux function as a coordinate. Next, we verify these results for the virtual-casing principle in standard cylindrical coordinates.

4.2. Analytic reduction to a line integral

We analytically simplify the integral in (4.1) by integrating over the toroidal angle. When we do so, we introduce the complete elliptic integrals of the first and second kind, which are defined as

(4.2a,b)\begin{equation} K(k^2) = \int_0^{{\rm \pi}/2} \frac{\mathrm{d} \theta}{\sqrt{1 - k^2 \sin^2\theta}} \quad \text{and} \quad E(k^2) = \int_0^{{\rm \pi}/2} \sqrt{1 - k^2 \sin^2\theta} \,\mathrm{d} \theta, \end{equation}

respectively. Our main analytical result gives a formula for the vector potential of the virtual-casing principle as a one-dimensional integral.

Theorem 4.1 Let $\varGamma$ be a surface of revolution with a generating curve $\gamma$ that satisfies the assumptions of § 2.1. Consider the vector potential

(4.3)\begin{equation} \boldsymbol{A}_S(\boldsymbol{x}) ={-}\frac 1{4{\rm \pi}} \iint_\varGamma \frac{\boldsymbol{n}(\boldsymbol{y}) \times \boldsymbol{B}^{\rm pol}(\boldsymbol{y})} {\left\| \boldsymbol{x} - \boldsymbol{y} \right\|} \,\mathrm{d} \varGamma(\boldsymbol{y}) \end{equation}

for $\boldsymbol {x} = (R,\phi _0,Z) = (r(t_0), \phi _0, z(t_0))$ in cylindrical coordinates. Define the quantities

(4.4a,b)\begin{equation} \alpha = \alpha(t; \boldsymbol{x}) = r(t)^2 + R^2 + (Z - z(t))^2 \quad \text{and} \quad \beta = \beta(t; \boldsymbol{x}) = 2 R r(t) \end{equation}

and set the modulus

(4.5)\begin{equation} k^2 = k(t; \boldsymbol{x})^2 = \frac{2\beta}{\alpha + \beta} = \frac{4 R r(t)}{(R + r(t))^2 + (Z - z(t))^2}. \end{equation}

Then the vector potential can be expressed as

(4.6)\begin{equation} \boldsymbol{A}_S(\boldsymbol{x}) ={-}\boldsymbol{e}_\phi(\phi_0) \int_0^L \mathcal{A}(t; \boldsymbol{x}) \,\mathrm{d} t \end{equation}

with the scalar-valued integrand

(4.7)\begin{equation} \mathcal{A}(t; \boldsymbol{x}) = \frac{1}{4{\rm \pi}} \left( \frac{4}{\sqrt{\alpha + \beta}} \right) \left[\dfrac{\partial \psi}{\partial z} r'(t) - \dfrac{\partial \psi}{\partial r} z'(t) \right] \left( \frac{2}{k^2} \left(K(k^2) - E(k^2) \right) - K(k^2) \right). \end{equation}

Proof. The unit outward normal vector of $\varGamma$ at $\boldsymbol {y}$ is given (by our assumption on the orientation of $\gamma$) by

(4.8)\begin{align} \boldsymbol{n}(\boldsymbol{y}(\phi,t)) & =\left. \frac{\partial{\boldsymbol{y}}}{\partial{\phi}}\times\frac{\partial{\boldsymbol{y}}}{\partial{t}} \right/ \left\| \frac{\partial{\boldsymbol{y}}}{\partial{\phi}}\times\frac{\partial{\boldsymbol{y}}}{\partial{t}} \right\| \nonumber\\ & = \frac{z'(t) \boldsymbol{e}_r(\phi) - r'(t) \boldsymbol{e}_z} {\sqrt{r'(t)^2 + z'(t)^2}}. \end{align}

From the representation (3.1) of $\boldsymbol {B}^{{\rm pol}}$ in axisymmetry, it follows that

(4.9)\begin{equation} (\boldsymbol{n}(\boldsymbol{y}(\phi,t)) \times \boldsymbol{B}^{{\rm pol}}(\boldsymbol{y}(\phi,t))) \mathcal{J}(t) = \left[\dfrac{\partial \psi}{\partial z} r'(t) - \dfrac{\partial \psi}{\partial r} z'(t) \right] \boldsymbol{e}_\phi(\phi). \end{equation}

Next, we consider the difference

(4.10)\begin{align} \boldsymbol{x} - \boldsymbol{y}(\phi, t) & = [R\boldsymbol{e}_r(\phi_0) + Z\boldsymbol{e}_z] - [r(t)\boldsymbol{e}_r(\phi)+ z(t)\boldsymbol{e}_z] \nonumber\\ & = ( R\cos\phi_0 - r(t)\cos\phi ) \boldsymbol{e}_x + ( R\sin\phi_0 - r(t)\sin\phi ) \boldsymbol{e}_y + (Z- z(t))\boldsymbol{e}_z, \end{align}

which we have expressed in both cylindrical and rectangular coordinates. Recall that the unit vector $\boldsymbol {e}_z$ is identical in both coordinate systems, and the remaining rectangular unit vectors $\boldsymbol {e}_x$ and $\boldsymbol {e}_y$ are related to standard cylindrical unit vectors by

(4.11a,b)\begin{equation} \boldsymbol{e}_r(\phi) = \cos\phi \, \boldsymbol{e}_x + \sin\phi \, \boldsymbol{e}_y \quad \text{and} \quad \boldsymbol{e}_\phi(\phi) ={-}\sin\phi \, \boldsymbol{e}_x + \cos\phi \, \boldsymbol{e}_y. \end{equation}

From the representation in rectangular coordinates, we use the trigonometric identity $\cos (\phi _0-\phi )=\cos \phi _0\cos \phi + \sin \phi _0\sin \phi$ and immediately obtain

(4.12)\begin{align} \left\| \boldsymbol{x} - \boldsymbol{y}(\theta, t) \right\|^2 = R^2 + r(t)^2 + (Z - z(t))^2 - 2 R r(t) \cos(\phi_0-\phi) = \alpha - \beta \cos(\phi_0 - \phi). \end{align}

We have now shown that the surface integral (4.1) for the vector potential is equivalent to the following double integral over a rectangle in the parameter domain:

(4.13)\begin{equation} \boldsymbol{A}_S(\boldsymbol{x}) ={-}\frac{1}{4{\rm \pi}} \int_0^L \left[\dfrac{\partial \psi}{\partial z} r'(t) - \dfrac{\partial \psi}{\partial r} z'(t) \right] \int_0^{2{\rm \pi}} \frac{\boldsymbol{e}_\phi(\phi)}{\sqrt{\alpha(t) - \beta(t) \cos(\phi_0 - \phi)}} \,\mathrm{d} \phi \,\mathrm{d} t. \end{equation}

We may analytically evaluate the inner integral using trigonometric identities and known integral formulae. We use the identities

(4.14)\begin{equation} \begin{cases} \sin(\phi+\phi_0) = \sin\phi\cos\phi_0 + \cos\phi\sin\phi_0\\ \cos(\phi+\phi_0) = \cos\phi\cos\phi_0 - \sin\phi\sin\phi_0 \end{cases} \end{equation}

to compute that

(4.15)\begin{align} \int_0^{2{\rm \pi}} \frac{-\sin\phi \, \boldsymbol{e}_x + \cos\phi \, \boldsymbol{e}_y} {\sqrt{\alpha - \beta\cos(\phi_0 - \phi)}} \,\mathrm{d} \phi & = \int_{0}^{2{\rm \pi}} \frac{-\sin(\phi+\phi_0) \, \boldsymbol{e}_x + \cos(\phi+\phi_0) \, \boldsymbol{e}_y} {\sqrt{\alpha - \beta\cos\phi}} \,\mathrm{d} \phi \nonumber\\ & = \boldsymbol{e}_\phi(\phi_0) \int_0^{2{\rm \pi}} \frac{\cos\phi \, \mathrm{d} \phi}{\sqrt{\alpha - \beta\cos\phi}} - \boldsymbol{e}_r(\phi_0) \int_0^{2{\rm \pi}} \frac{\sin\phi \, \mathrm{d} \phi}{\sqrt{\alpha - \beta\cos\phi}} \nonumber\\ & = \boldsymbol{e}_\phi(\phi_0) \int_0^{2{\rm \pi}} \frac{\cos\phi \, \mathrm{d} \phi}{\sqrt{\alpha - \beta\cos\phi}}. \end{align}

Here, we have used the fact that

(4.16)\begin{equation} \int_0^{2{\rm \pi}} \frac{\sin\phi \, \mathrm{d} \phi}{\sqrt{\alpha - \beta \cos\phi}} = 0, \end{equation}

which holds because it is the integral of an odd function over a single period.

Now, the remaining integral can be expressed in terms of the complete elliptic integrals, through the identity

(4.17)\begin{align} \int_0^{2{\rm \pi}} \frac{\cos\phi \, \mathrm{d} \phi}{\sqrt{\alpha - \beta\cos\phi}} & = 2\int_{-{\rm \pi}/2}^{{\rm \pi}/2} \frac{\cos (2\phi+{\rm \pi}) \, \mathrm{d} \phi}{\sqrt{\alpha - \beta\cos (2\phi+{\rm \pi})}} \nonumber\\ & = 2\int_{-{\rm \pi}/2}^{{\rm \pi}/2} \frac{-(1 - 2\sin^2\phi) \, \mathrm{d} \phi}{\sqrt{\alpha + \beta(1 - 2\sin^2\phi)}} \nonumber\\ & = \frac{-4}{\sqrt{\alpha+\beta}} \int_0^{{\rm \pi}/2} \frac{1-2\sin^2\phi} {\sqrt{1 - k^2 \sin^2\phi}} \mathrm{d} \phi \nonumber\\ & = \frac{-4}{\sqrt{\alpha+\beta}} \left( K(k^2) - \frac 2{k^2} (K(k^2) - E(k^2))\right). \end{align}

We obtain the last equality (4.17) from the definition of $K$ and from a formula of Gradshteyn & Ryzhik (2014, 2.584-4). The desired result now immediately follows.

The Kapur–Rokhlin quadrature rule that we analysed in § 2.3.2 applies to integrands with a logarithmic singularity. With the result we just obtained, we can verify, by analysing the behaviour of the complete elliptic integrals, that the integrand $\mathcal {A}$ is indeed logarithmically singular as $t \to t_0$, in agreement with the results from Chance (1997, § V). Specifically, as the modulus $k$ tends to $1$, the second-kind integral $E$ is continuous and bounded, and the first-kind elliptic integral $K$ is logarithmically singular. The mapping $t \mapsto k(t; \boldsymbol {x})^2$ is continuous, so $E(k(t;\boldsymbol {x})^2)$ is continuous and $K(k(t;\boldsymbol {x})^2)$ is logarithmically singular as $t \to t_0$. Readers interested in more detail regarding these results are referred to Lemma A.2 in the Appendix.

5. A Kapur–Rokhlin scheme for the virtual-casing principle

We may compute the vector potential by (4.6) of Theorem 4.1 using the Kapur-Rokhlin quadrature rule for periodic functions introduced in § 2.3.2. Let $\boldsymbol {x} \in \varGamma$ be given, corresponding to parameters $(\phi _0,t_0)$. We reiterate that a univariate integral expression for the vector potential is

(5.1)\begin{equation} \boldsymbol{A}_S(\boldsymbol{x}) ={-}\boldsymbol{e}_\phi(\phi_0) \int_0^L \mathcal{A}(t; \boldsymbol{x}) \,\mathrm{d} t. \end{equation}

The Kapur–Rokhlin quadrature of § 2.3.2 applies directly because the integration interval $[0,L]$ is identical, by periodicity, to the symmetric interval $[t_0-L/2, t_0+L/2]$ and because the integrand $\mathcal {A}(t; \boldsymbol {x})$ is logarithmically singular as $t \to t_0$.

Given $M \in \mathbb {N}$, we generate $2M-1$ quadrature nodes $t_i = t_0 + ih$ for $i = \pm 1, \dots, \pm (M-1), M$ with spacing $h=L/(2M)$. We evaluate $\mathcal {A}_i = \mathcal {A}(t_i; \boldsymbol {x})$, and it follows that

(5.2)\begin{equation} \int_0^L \mathcal{A}(t; \boldsymbol{x}) \,\mathrm{d} t = h \left[ \sum_{1 \leqslant |i| \leqslant M-1} \mathcal{A}_i + \mathcal{A}_M \right] + h \sum_{1 \leqslant |j| \leqslant n} \gamma_j (\mathcal{A}_j + \mathcal{A}_{{-}j}) + O(h^n). \end{equation}

As periodicity has removed our considerations about expanding the integration domain, the parameter $n$ can be taken as large as we like (as long as $M \geqslant n$ to define the quadrature rule). However, in practice this Kapur–Rokhlin scheme is known to be unstable for $n$ larger than about 10 because the weights $\gamma _j$ are sign-indefinite and grow large in magnitude.

6. Numerical results

In this section, we illustrate the Kapur–Rokhlin quadrature schemes from §§ 2.3.2 and 5 in two different calculations. Throughout, we test the schemes for an axisymmetric plasma boundary given by the level set $\psi =0$ of the poloidal flux function given by

(6.1)\begin{equation} \psi(r,z) = \frac{\kappa F_B}{2 R_0^3 q_0} \left[\frac 14 (r^2 - R_0^2)^2 + \frac 1{\kappa^2} r^2 z^2 - a^2 R_0^2 \right]. \end{equation}

This flux function is a solution to the Grad–Shafranov equation with the Solov'ev profiles $\mu _{0}p(\psi )=-[F_{B}(\kappa + 1/\kappa )/( R_0^{3}q_0)]\psi$ and $F(\psi )=F_B$, where $p(\psi )$ is the plasma pressure profile, and $F(\psi )=rB_{\phi }$, with $B_\phi$ the toroidal magnetic field (Lütjens, Bondeson & Sauter 1996; Lee & Cerfon 2015). The parameters $R_0$ and $q_0$ may be interpreted as the major radius and safety factor at the magnetic axis, and $\kappa$ and $a$ as the elongation and minor radius of the plasma boundary. All numerical tests in this article use the fusion relevant values $F_B = R_0 = q_0 = 1$ and $\kappa = 1.7$ and $a = 1/3$. The level set $\psi = 0$ may be parameterised by the functions (Lütjens et al. 1996; Lee & Cerfon 2015)

(6.2a,b)\begin{equation} (r(t))^2 = R_0^2 + 2aR_0 \cos t \quad \text{and} \quad z(t) = \kappa a \frac{R_0}{r(t)}\sin t \end{equation}

for $t \in [0, 2{\rm \pi} ]$.

6.1. Double-layer identity

For our first numerical verification, we consider an identity associated with harmonic functions. Consider the Green's function $G(\boldsymbol {x}, \boldsymbol {y}) = (4 {\rm \pi}\left \| \boldsymbol {x} - \boldsymbol {y} \right \|)^{-1}$ for Laplace's equation in three dimensions. It satisfies the double-layer jump condition (Malhotra et al. 2019b)

(6.3)\begin{equation} \iint_\varGamma \frac{\partial G(\boldsymbol{x}, \boldsymbol{y})}{\partial \boldsymbol{n}(\boldsymbol{y})} \,\mathrm{d} \varGamma(\boldsymbol{y}) = \frac 1{4{\rm \pi}} \iint_\varGamma \frac{\boldsymbol{n}(\boldsymbol{y}) \boldsymbol{ \cdot} (\boldsymbol{x} - \boldsymbol{y})} {\left\| \boldsymbol{x} - \boldsymbol{y} \right\|^3} \,\mathrm{d} \varGamma(\boldsymbol{y}) ={-} \frac 12 \end{equation}

for $\boldsymbol {x} \in \varGamma$. It follows that

(6.4)\begin{equation} 1 + \frac 1{2{\rm \pi}} \iint_\varGamma \frac{\boldsymbol{n}(\boldsymbol{y}) \boldsymbol{ \cdot} (\boldsymbol{x} - \boldsymbol{y})} {\left\| \boldsymbol{x} - \boldsymbol{y} \right\|^3} \,\mathrm{d} \varGamma(\boldsymbol{y}) = 0, \end{equation}

again for $\boldsymbol {x} \in \varGamma$. Following identical methodology to the proof of Theorem 4.1, we integrate out the toroidal angle analytically to obtain the one-dimensional integral identity

(6.5)\begin{align} 1 + \frac 1{2{\rm \pi}} \int_0^L \frac{4r}{(\alpha+\beta)^{3/2}} \left\{ -\frac{2 z' R}{k^2} K(k^2) + \left( \frac{2 z' R}{k^2} + \frac{z'(R-r) - r'(Z-z)}{1-k^2} \right) E(k^2) \right\} \mathrm{d} t = 0. \end{align}

In this expression, we have suppressed the dependence of $\{r,z,r',z',\alpha,\beta,k^2\}$ on $t$ for clarity. As usual, we have also used the identification $\boldsymbol {x} = (R, \phi _0, Z) = (r(t_0), \phi _0, z(t_0))$. The integrand is logarithmically singular because of the presence of the singular elliptic integral $K(k^2)$, and because the other seemingly singular coefficient is actually bounded when $r''(t_0)$ and $z''(t_0)$ exist (which is the case in our example), with

(6.6)\begin{equation} \lim_{t \to t_0} \frac{z'(t)(R-r(t)) - r'(t)(Z-z(t))}{1-k(t)^2} = 2 R^2 \left( \frac{r''(t_0) z'(t_0) - r'(t_0) z''(t_0)}{r'(t_0)^2 + z'(t_0)^2} \right). \end{equation}

We conclude that the integral in (6.5) is of the Kapur–Rokhlin form from § 2.3.2.

Figure 3 illustrates the performance of the 10th-order periodic Kapur–Rokhlin quadrature scheme for verifying the identity (6.5). We compare this method with the alternating trapezoidal rule, which uses common quadrature weight $h={\rm \pi} /M$ and quadrature nodes $\tilde {t}_i = t_0 + (i - \frac 12) h$ for $i=0, \pm 1, \dots, \pm (M-1), M$ that straddle the singularity at $t_0$. We find that the Kapur–Rokhlin scheme achieves the theoretical 10th-order accuracy, and that we obtain the full accuracy of the quadrature scheme using about 175 quadrature nodes. Because the behaviour of the integrand is not generally symmetric about $t=t_0$, the alternating trapezoidal rule performs poorly. In order to avoid unfairly representing the performance of singularity subtraction schemes currently used in plasma physics applications, we have not attempted to code our own versions of them to compare their accuracy and their convergence rate in this figure. However, we explain in Malhotra et al. (2019b) why all these singularity subtraction schemes are expected have low-order convergence, with second-order convergence expected for most of them.

Figure 3. Error convergence for the double-layer identity (6.5) at $t_0 = 1$.

6.2. Virtual-casing principle

For our second test, we evaluate the accuracy of our method for calculating the normal component of the magnetic field due to the plasma current on the plasma boundary, i.e. $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {B}^{{\rm pol}}_V$ on $\varGamma$. The plasma equilibrium we consider is the Solov'ev equilibrium described previously (Lütjens et al. 1996; Lee & Cerfon 2015). As we do not know the analytic solution to this problem, we compare our implementation with an existing high-order accurate implementation from Malhotra et al. (2019b). The method used in that implementation is different from the approach presented here in several ways, making it appropriate for our verification. Specifically, the code presented in Malhotra et al. (2019b) views the plasma equilibrium as a fully three-dimensional equilibrium, and does not assume axisymmetry. Furthermore, Malhotra et al. (2019b) obtain $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {B}^{{\rm pol}}_V$ on $\varGamma$ directly via a direct evaluation of (3.3), as opposed to first computing the vector potential. Finally, the Cauchy principal value of (3.3) is numerically evaluated via a partition of unity scheme to handle the singulary of the integrand. We take the result from a high-resolution calculation of Malhotra et al. (2019b) for this problem as the ground truth, against which we test the accuracy of our approach as a function of the number of quadrature points.

Figure 4 illustrates our results when using the 10th-order Kapur–Rokhlin scheme for this problem. Again, we find that the Kapur–Rokhlin scheme achieves the theoretical convergence rate, and that the accuracy has converged by about 400 quadrature nodes. In both this example and the double-layer identity example, we observe that the Kapur–Rokhlin scheme errors do not converge all the way to machine precision. This is one known drawback of the Kapur–Rokhlin methods, driven partially by the instabilities caused by the correction weights. Nonetheless, in contexts where full machine precision is not necessary, this scheme provides high accuracy and is easy to implement by reading the correction weights $\gamma _j$ from a table of pre-computed values.

Figure 4. Comparison with existing code to compute the normal component of the magnetic field.

7. Conclusion

For axisymmetric confinement fusion systems, a direct implementation of the Kapur–Rokhlin quadrature scheme yields high accuracy for the evaluation of the layer potentials commonly encountered in magnetostatic and magnetohydrodynamic calculations. Using the table of quadrature weights given in the original article by Kapur & Rokhlin (1997), the scheme is as easy to implement as the trapezoidal rule, and, unlike commonly used methods, does not require any manipulation of the singular integrands. We demonstrated how to implement it for the evaluation of a double-layer potential and for the virtual-casing principle, obtaining 10th-order convergence in both cases.

At the moment, our approach is restricted to singular integrals along smooth boundaries, and cannot be applied to magnetic surfaces with one or several X-points. The development of an efficient and accurate quadrature scheme for this importance case is the subject of ongoing work, with results to be reported at a later date.