Adaptive Hybridizable Discontinuous Galerkin discretization of the Grad-Shafranov equation by extension from polygonal subdomains

TL;DR

This paper introduces an adaptive hybridizable discontinuous Galerkin method for solving the Grad-Shafranov equation on curved domains, using polygonal subdomains and automatic mesh refinement to accurately model plasma equilibria.

Contribution

It presents a novel high-order solver that avoids complex geometry-conforming meshes by employing a transfer technique on polygonal subdomains with adaptive mesh refinement.

Findings

Effective error estimator for plasma equilibria

Automatic boundary updating maintains geometric accuracy

Suitable for realistic plasma confinement geometries

Abstract

We propose a high-order adaptive numerical solver for the semilinear elliptic boundary value problem modelling magnetic plasma equilibrium in axisymmetric confinement devices. In the fixed boundary case, the equation is posed on curved domains with piecewise smooth curved boundaries that may present corners. The solution method we present is based on the hybridizable discontinuous Galerkin method and sidesteps the need for geometry-conforming triangulations thanks to a transfer technique that allows to approximate the solution using only a polygonal subset as computational domain. Moreover, the solver features automatic mesh refinement driven by a residual-based a posteriori error estimator. As the mesh is locally refined, the computational domain is automatically updated in order to always maintain the distance between the actual boundary and the computational boundary of the order of the local mesh diameter. Numerical evidence is presented of the suitability of the estimator as an approximate error measure for physically relevant equilibria with pressure pedestals, internal transport barriers, and current holes on realistic geometries.

Tables

Table 1. Table 1: Estimated convergence rates (e.c.r.) for the approximate solutions ψ h subscript 𝜓 ℎ \psi_{h} and 𝒒 h subscript 𝒒 ℎ \boldsymbol{q}_{h} as well as for the error estimator η 𝜂 \eta and all its separate components in the case of uniform refinement. The exact Solov’ev solution is given in Equation ( 22 ). For a given quantity U 𝑈 U we denote the L 2 superscript 𝐿 2 L^{2} error at the k − limit-from 𝑘 k- th refinement level by E ​ ( U ) 2 k 𝐸 subscript superscript 𝑈 𝑘 2 E(U)^{k}_{2} . Then, the estimated convergence rate is computed through the formula e.c.r = log 2 ⁡ ( E ​ ( U ) 2 k / E ​ ( U ) 2 k + 1 ) e.c.r subscript 2 𝐸 subscript superscript 𝑈 𝑘 2 𝐸 subscript superscript 𝑈 𝑘 1 2 \text{e.c.r}=\log_{2}{\left(E(U)^{k}_{2}/E(U)^{k+1}_{2}\right)} .

	Polynomial degree $k=1$								Polynomial degree $k=2$
	${\psi }_h$	${𝒒}_h$	$\eta$	${\eta }_1$	${\eta }_2$	${\eta }_3$	${\eta }_4$	${\eta }_5$	${\psi }_h$	${𝒒}_h$	$\eta$	${\eta }_1$	${\eta }_2$	${\eta }_3$	${\eta }_4$	${\eta }_5$
$h\to h/2$	1.32	0.79	1.66	1.68	1.64	1.49	2.1	1.81	2.57	1.66	2.01	2.04	2.02	1.63	1.52	2.78
$h/2\to h/4$	1.95	2.34	1.8	1.84	1.62	1.59	2.05	1.5	2.7	2.89	2.63	2.64	2.81	2.55	2.93	3.06
$h/4\to h/8$	1.97	2.11	1.87	1.89	1.82	1.78	2.28	1.81	2.8	2.26	2.66	2.66	2.75	2.62	2.56	2.81
$h/8\to h/16$	1.99	2.13	1.93	1.95	1.93	1.86	2.33	1.94	2.97	3.02	2.88	2.88	2.91	2.83	2.74	2.85
	Polynomial degree $k=3$								Polynomial degree $k=4$
	${\psi }_h$	${𝒒}_h$	$\eta$	${\eta }_1$	${\eta }_2$	${\eta }_3$	${\eta }_4$	${\eta }_5$	${\psi }_h$	${𝒒}_h$	$\eta$	${\eta }_1$	${\eta }_2$	${\eta }_3$	${\eta }_4$	${\eta }_5$
$h\to h/2$	3.46	1.66	1.64	1.61	1.82	1.99	1.96	2.82	3.94	2.14	1.39	1.4	1.65	1.28	1.64	2.47
$h/2\to h/4$	3.52	3.62	3.53	3.53	3.58	3.53	3.91	3.93	4.37	4.02	4.34	4.34	4.29	4.21	4.53	4.56
$h/4\to h/8$	3.13	2.3	3.26	3.26	3.2	3.34	3.44	3.48	3.64	2.7	3.88	3.89	3.63	3.81	3.92	3.96
$h/8\to h/16$	4.26	3.78	3.83	3.83	3.81	3.85	3.99	3.71	5.21	4.56	4.78	4.79	4.68	4.67	4.79	4.79

Equations

$$-{\Delta }^{\ast }\psi ={\mu }_0{r}^2\frac{dp}{d\psi }+\frac{1}{2}\frac{d{g}^2}{d\psi }=:F(r,\psi ).$$

$${\Delta }^{\ast }\psi :=r^{2}div\left(\frac{1}{r^2}grad\psi \right)=r{\partial }_r\left(\frac{1}{r}{\partial }_r\psi \right)+{\partial }_z^2\psi =r{\partial }_r\left(\frac{1}{r}{\partial }_r\psi \right)+r{\partial }_z\left(\frac{1}{r}{\partial }_z\psi \right)=r\tilde{\nabla }\cdot \left(\frac{1}{r}\tilde{\nabla }\psi \right).$$

$$\mathbf{B}=\frac{1}{r}\nabla \psi \times {\mathbf{e}}_{\varphi }+\frac{g(\psi )}{r}{\mathbf{e}}_{\varphi }$$

$${\mu }_0\mathbf{J}=\frac{1}{r}\frac{dF}{d\psi }\nabla \psi \times {\mathbf{e}}_{\varphi }-\frac{1}{r}{\Delta }^{\ast }\psi {\mathbf{e}}_{\varphi }$$

$$-r\tilde{\nabla }\cdot \left(\frac{1}{r}\tilde{\nabla }\psi \right)=$$

$$\psi =$$

$$\mathbf{q}-\frac{1}{r}\tilde{\nabla }\psi =$$

$$-\tilde{\nabla }\cdot \mathbf{q}=$$

$$\psi =$$

$${\mathsf{T}}^h:=\left\{K\in {\mathcal{T}}^h:\overline{K}\subset \Omega \right\},\partial {\mathsf{T}}^h:=\left\{\partial K:K\in {\mathsf{T}}^h\right\},\text{and}{\Omega }^h:={\left({\cup }_{K\in {\mathsf{T}}^h}\overline{K}\right)}^{\circ }.$$

$${\mathcal{E}}^{\partial }:=\left\{e\in {\mathcal{E}}^h:e\subset {\Gamma }^h\right\}\text{ and }{\mathcal{E}}_0^{\circ }:=\left\{e\in {\mathcal{E}}^h:e\text{nsubset}{\Gamma }^h\right\}$$

$${\mathsf{T}}_c^h:=\left\{K\in {\mathcal{T}}^h|K\cap \overline{\Omega }\ne \varnothing \right\},$$

$$\psi (\mathbf{x})=\psi (\bar{\mathbf{x}})-{\int }_0^{d(\mathbf{x},\overline{\mathbf{x}})}r(\mathbf{x}+\mathbf{t}(\mathbf{x})s)\mathbf{q}(\mathbf{x}+\mathbf{t}(\mathbf{x})s)\cdot \mathbf{t}(\mathbf{x})ds,$$

$$E(p):K^e\cup K_{ext}^e$$

$$\mathbf{x}$$

$${\phi }_h(\mathbf{x}):=\psi (\bar{\mathbf{x}})-{\int }_0^{d(\mathbf{x},\overline{\mathbf{x}})}r(\mathbf{x}+\mathbf{t}(\mathbf{x})s)E({\mathbf{q}}_h)(\mathbf{x}+\mathbf{t}(\mathbf{x})s)\cdot \mathbf{t}(\mathbf{x})ds$$

$$\mathbf{q}-\frac{1}{r}\tilde{\nabla }\psi =$$

$$-\tilde{\nabla }\cdot \mathbf{q}=$$

$$\psi =$$

$$[[\mathbf{q}]]=$$

$$\psi =$$

$$[[\mathbf{q}]]:={\mathbf{q}}^{+}\cdot {\mathbf{n}}^{+}+{\mathbf{q}}^{-}\cdot {\mathbf{n}}^{-}$$

$$\hat{\mathbf{q}}\cdot \mathbf{n}:=\mathbf{q}\cdot \mathbf{n}-\tau (\psi -\hat{\psi }).$$

$${\mathbf{V}}_h$$

$$W_h$$

$$M_h$$

$$(\cdot ,\cdot {)}_{{\mathsf{T}}^h}:=\sum _{K\in {\mathsf{T}}^h}(\cdot ,\cdot {)}_K⟨\cdot ,\cdot {⟩}_{\partial {\mathsf{T}}^h}:=\sum _{K\in {\mathsf{T}}^h}⟨\cdot ,\cdot {⟩}_{\partial K},$$

$$(r{\mathbf{q}}_h,\mathbf{v}{)}_{{\mathsf{T}}^h}+({\psi }_h,\nabla \cdot \mathbf{v}{)}_{{\mathsf{T}}^h}$$

$$-(\nabla \cdot {\mathbf{q}}_h,w{)}_{{\mathsf{T}}^h}-⟨\tau {\psi }_h,w{⟩}_{\partial {\mathsf{T}}^h}$$

$$⟨{\mathbf{q}}_h\cdot \mathbf{n},\mu {⟩}_{\partial {\text{T}}^h\setminus {\Gamma }^h}+⟨\tau {\psi }_h,\mu {⟩}_{\partial {\text{T}}^h\setminus {\Gamma }^h}$$

$$⟨{\phi }_h,\mu {⟩}_{{\Gamma }^h}$$

$${\phi }_h(\mathbf{x})=-{\int }_0^{d(\mathbf{x},\overline{\mathbf{x}})}r(\mathbf{x}+\mathbf{t}(\mathbf{x})s)E({\mathbf{q}}_h)(\mathbf{x}+\mathbf{t}(\mathbf{x})s)\cdot \mathbf{t}(\mathbf{x})ds.$$

$$[r{]}_{ij}$$

$$[\tau {]}_{ij}$$

$$\displaystyle\left[\begin{array}[]{cc}r&\nabla\cdot^{\top}\\ -\,\nabla\cdot&-\tau\end{array}\right]\,\left[\begin{array}[]{c}\boldsymbol{q}_{h}\\ \psi_{h}\end{array}\right]$$

$$\displaystyle\left[\begin{array}[]{c}0\\ \varphi_{h}\end{array}\right]+\left[\begin{array}[]{cc}\boldsymbol{n}\,\mathbf{I}^{\circ}_{e}&-\tau\,\mathbf{I}^{\circ}_{e}\\ 0&0\end{array}\right]\,\left[\begin{array}[]{c}\boldsymbol{q}_{h}\\ \psi_{h}\end{array}\right]$$

$$\left[\begin{array}[]{c}0\\ \varphi_{h}\end{array}\right]=\left[\begin{array}[]{cc}0&0\\ -\mathbf{Q}&0\end{array}\right]\left[\begin{array}[]{c}\boldsymbol{q}_{h}\\ \psi_{h}\end{array}\right]$$

$$[\tilde{\mathbf{Q}}{]}_{ij}:=⟨{\int }_0^{d(\mathbf{x},\overline{\mathbf{x}})}r(\mathbf{x}+\mathbf{t}(\mathbf{x})s)E({\mathbf{\varphi }}_i)(\mathbf{x}+\mathbf{t}(\mathbf{x})s)\cdot \mathbf{t}(\mathbf{x})ds,{\mu }_i(\mathbf{x}){⟩}_{\partial K}.$$

Full text

fn1fn1footnotetext: A.J.C. and T. S-V have been partially funded by the US Department of Energy. Grant No. DE-FG02-86ER53233.fn2fn2footnotetext: M. S. has been partially funded by CONICYT–Chile through FONDECYT project No. 1160320 and by Project AFB170001 of the PIA Program: Concurso Apoyo a Centros Cientificos y Tecnologicos de Excelencia con Financiamiento Basal.

1 Introduction

In toroidally axisymmetric configurations and in the absence of flows, the steady-state equations of magnetohydrodynamics (MHD) yield the following partial differential equation for the poloidal flux function $\psi=\psi(r,z)$, known as the Grad-Shafranov equation [1, 2, 3]

[TABLE]

In Equation (1), $r$ is the radial coordinate in the $(r,\phi,z)$ coordinate system naturally associated with the toroidal geometry, $\phi$ being the ignorable coordinate for the axisymmetric configurations we consider here, $\psi=\psi_{p}/2\pi$ where $\psi_{p}$ is the poloidal magnetic flux, $p=p(\psi)$ is the plasma pressure, $2\pi g(\psi)=-I_{p}$ is the net poloidal current flowing in the plasma and the toroidal field coils, and the elliptic toroidal operator is defined by

[TABLE]

Above, for simplicity in the manipulations we have defined the operator $\widetilde{\nabla}:=(\partial_{r},\partial_{z})$ that acts formally like a vector of partial derivatives independent of the coordinate system.

Both $p$ and $g$ are free functions of $\psi$, which are determined from other physical processes or experimental data, and taken as input to the partial differential equation. In general, $p$ and $g$ are such that $F(r,\psi)$ is a nonlinear function of $\psi$, so that the Grad-Shafranov equation is a semi-linear partial differential equation. Together with the boundary conditions, $p$ and $g$ determine the nature of the MHD equilibrium. Once $\psi$ is computed, the equilibrium magnetic configuration is fully determined, through the relations

[TABLE]

where $\mathbf{B}$ is the magnetic field and $\mathbf{J}$ is the current density, and $\mathbf{e}_{\phi}$ is the unit vector in the toroidal direction.

In this article, we will focus on fixed boundary equilibria, for which the boundary $\Gamma$ of the computational domain $\Omega$ is known, and corresponds to the boundary of the confinement region of the plasma. Physically, it must be a level set of $\psi$, and without loss of generality, we can let $\Gamma:=\partial\Omega$ be the level set $\psi=0$. We are therefore interested in the following Dirichlet boundary value problem

[TABLE]

For plasma physics applications, the solution $\psi$ to (4), the corresponding magnetic field $\mathbf{B}$ and the current density $\mathbf{J}$ given by Equations (2) and (3), are used as input for stability, transport, and radio-frequency (RF) wave propagation and heating solvers [4, 5, 6, 7, 8, 9, 10]. This gives stringent performance requirements on numerical solvers for the Grad-Shafranov equation. The solver should be fast, so that the time to compute the equilibrium configuration is negligible compared to the run time of the stability, transport, or RF wave solvers. The solver should also be accurate, because some of the physical quantities of interest depend not only on $\psi$, but also on the first derivatives of $\psi$, as is obvious for $\mathbf{B}$ and $\mathbf{J}$, and on the second derivatives of $\psi$, as is the case for the magnetic curvature for example. Because of the central role of the Grad-Shafranov equation in magnetic confinement fusion, many numerical solvers for Eq.(4) have been developed in the last decades, relying on a vast range of formulations and numerical schemes. A good summary of early efforts can be found in [11]. More recently, approaches based on bi-cubic finite elements [7, 12, 13], on spectral elements [6, 14], on the hybridizable discontinuous Galerkin method[15], and on integral equation methods [16, 17] have led to fast, high order accurate, and flexible solvers. Even so, none of these solvers simultaneously satisfy the four criteria required for optimal performance in magnetic confinement fusion applications, which can be listed as follows: 1) the solver must be fast; 2) it must be able to handle arbitrary boundaries $\Gamma$, which may or may not have corners (corresponding to magnetic X-points); 3) it must compute derivatives with high accuracy; and 4) it must have automatic adaptive refinement capabilities, in order to resolve the strong gradients in internal transport barriers and in edge pedestals without having a fine grid throughout the computational domain, where the solution is typically very smooth and requires few grid points for good accuracy. In this article, we present the first numerical Grad-Shafranov solver which satisfies these four performance requirements.

To achieve this goal, we added adaptive refinement capabilities to the Hybridizable Discontinuous Galerkin (HDG) Grad-Shafranov solver we originally presented in [15]. The solver relies on reformulating the problem as a first order system [18], which takes the form

[TABLE]

The auxiliary variable $\boldsymbol{q}$ will be referred to as the flux. This mixed formulation has the practical advantage of discretizing $\boldsymbol{q}$ directly, thus providing additional accuracy for the the physically meaningful quantity. Our enhanced solver exploits the natural suitability of discontinuous Galerkin methods for parallel computation, sidesteps the geometrical complexities by carrying on the computations on a polygonal subdomain of $\Omega$ discretized by a simple non-fitting and shape regular triangulation, handling the curved boundaries with a high order transferring technique, approximates the partial derivatives of $\psi$ directly, and features automatic mesh refinement driven by a local error estimator.

The structure of the paper is as follows. Starting with geometric considerations, Section 2 describes the discretization of the mixed form (5), introducing the basic concepts of the hybridizable discontinuous Galerkin method and the components of the solution process including the treatment of curved boundaries, the accelerated iteration process to handle the non-linearity of the Grad-Shafranov equation, and a post-processing step yielding an approximation to $\psi_{h}$ that converges with an additional order of accuracy. Throughout Section 2 the problem is posed on a fixed computational domain with a given mesh. Adaptivity is then addressed in Section 3, which starts by introducing a residual-based local error estimator and by briefly discussing choices for the element marking strategy. The remaining part of the section discusses the issue of refining the embedded triangulation while maintaining the distance between the computational domain and the boundary always on the order of the local mesh parameter. The resulting strategy generates a sequence of updates of the computational domain that approximate the physical domain by exhaustion as the refinement progresses. In Section 4 we present numerical experiments to demonstrate the efficiency and reliability of the error estimator, as well as the convergence properties of the numerical solution. The experiments are carried out in realistic geometries first for a Solov’ev equilibrium for which an exact solution is available, and then for physically relevant benchmark problems with sharp and localized features, specifically an equilibrium with a pressure pedestal, and equilibrium with an internal transport barrier, and an equilibrium with a deep current hole. Concluding remarks are given in the final Section 5.

2 The discrete problem

Before describing the adaptive algorithm, we will discuss the problem on a fixed, uniform and embedded polygonal mesh. Most of the details have been described in [15]. However, starting with this standard case will allow us to introduce the notation and the fundamental ideas underpinning our HDG approach, as well as our treatment of curved boundaries and the iterative method to treat the non-linearity of the equation. Moreover, the solution of the problem in this setting will constitute the starting point for the adaptive algorithm. It is therefore worth repeating the key elements of the numerical scheme for the standard situation here.

The HDG formulation of the problem depends on the specific spatial discretization of the domain where the equation is posed. Therefore, we will first describe the choice of polygonal subdomain where the problem will be discretized. The use of a polygonal subset of $\Omega$ as the computational domain creates the need to communicate the Dirichlet boundary conditions from the “true” boundary to that of the polygonal subdomain where the computations are carried out. Hence, we will then describe the high order transfer scheme that will be used to impose the boundary conditions on the computational domain. With all these ingredients in place, it will be then possible to pose the discrete problem and describe the numerical method in detail. A rigorous analysis of the method for general semilinear elliptic equations is the subject of ongoing work [19].

2.1 The computational domain

We will start by defining the polygonal subdomain where the discrete system will be posed, henceforward the computational domain, and the grids that will be used for approximation. Following [20], the computational domain $\Omega^{h}$ will be chosen to be a polygonal subdomain of $\Omega$ obtained from a regular background triangulation as follows.

Consider $\mathcal{T}^{h}$ to be a triangulation of a polygonal domain containing $\overline{\Omega}$ and consisting of uniformly shape-regular triangles $K$ as in Figure 1 (left). The computational mesh $\mathsf{T}^{h}$ and the computational domain $\Omega^{h}$ (Figure 1 center) are, respectively, the set of triangles completely contained in $\Omega$ and its interior. More precisely, we define

[TABLE]

The boundary of the computational domain will be denoted by $\Gamma^{h}$ and set of all edges $e$ commonly referred to as the skeleton of the triangulation will be denoted by $\mathcal{E}^{h}$. We note that the skeleton can be decomposed as $\mathcal{E}^{h}=\mathcal{E}^{\partial}\cup\mathcal{E}^{\circ}$ where

[TABLE]

are the set of boundary edges and interior edges respectively. In addition, a companion grid consisting of those elements in $\mathcal{T}^{h}$ that constitute a minimal cover of $\overline{\Omega}$ will be defined

[TABLE]

the companion mesh consists of all the elements in $\mathsf{T}^{h}$ together with those background elements that intersect with the boundary of $\Omega$, as depicted on the right end of Figure 1. The need for this additional companion mesh will become apparent in Section 3.3, where the mesh refinement strategy will be discussed.

2.2 Extension from subdomains

The definition of the computational domain $\Omega^{h}$ as an unfitted and embedded subdomain seems to leave open the problem of defining the approximate solutions in the intermediate region $\Omega\setminus\overline{\Omega^{h}}$ corresponding to the area between the “true” boundary and the computational boundary. In particular, one must deal with the problem of defining the boundary conditions on $\Gamma^{h}$. This in fact can be dealt with in a natural way through the following transfer procedure.

Consider a point $\overline{\boldsymbol{x}}=(\overline{r},\overline{z})$ on the true boundary $\Gamma$ and a point $\boldsymbol{x}=(r,z)$ on the computational boundary $\Gamma^{h}$. We will denote by $\boldsymbol{t}$ the normalized vector anchored at $\boldsymbol{x}$ pointing towards $\overline{\boldsymbol{x}}$ and by $\sigma_{\boldsymbol{t}}(\boldsymbol{x})$ the line segment connecting them, henceforward the transfer path – as in Figure 2 (right). Then, integrating equation (5a) along a transfer path with direction vector given by $\boldsymbol{t}(\boldsymbol{x})$ it follows that

[TABLE]

where $d(\boldsymbol{x},\overline{\boldsymbol{x}})$ is the Euclidean distance between $\boldsymbol{x}$ and $\overline{\boldsymbol{x}}$. Therefore, as long as $\boldsymbol{q}$ is known along the transfer path, it is possible to represent exactly the value of $\psi$ at any point $\boldsymbol{x}$ of the computational boundary in terms of its value at one point $\overline{\boldsymbol{x}}$ of the physical boundary and the values of the flux. However, the value of $\boldsymbol{q}$ will be determined only within the computational domain $\Omega^{h}$ and thus we will resort to an approximation of $\boldsymbol{q}$ by extrapolation. In order to detail the extrapolation procedure we must first introduce the following extension of the computational domain.

Consider the set $\mathcal{Y}$ consisting of all endpoints of the edges $e\in\mathcal{E}^{\partial}_{0}$ and denote by $h_{loc}(\boldsymbol{y})$ the minimum diameter $h$ over all triangles containing $\boldsymbol{y}$. To every $\boldsymbol{y}\in\mathcal{Y}$ we will assign a unique point $\overline{\boldsymbol{y}}\in\Gamma$ and will denote by $\sigma_{t}(\boldsymbol{y})$ the straight line segment connecting them. The assignment must be done such that the Euclidean distance between them $d(\boldsymbol{y},\overline{\boldsymbol{y}})=\mathcal{O}(h_{loc}(\boldsymbol{y}))$, and that no two pair of such paths intersects. This can be done in different ways, for instance, following the procedure described in [20].

For a boundary edge $e\in\mathcal{E}^{\partial}_{0}$ we will denote by $K^{e}_{ext}$ the region enclosed by $\Gamma$, the paths corresponding to each endpoint of $e$, and the edge itself—as depicted on the right end of Figure 2. The union of all these patches covers the “un-meshed” gap $\Omega\setminus\overline{\Omega^{h}}$, as can be seen in the center of Figure 2. Each of the patches $K^{e}_{ext}$ can be unambiguously identified with the unique element $K^{e}$ with which they share the edge $e$. This allows us to define the extension $E(p):K^{e}\cup K^{e}_{ext}\rightarrow\mathbb{R}$ of a polynomial function $p:K^{e}\rightarrow\mathbb{R}$ as

[TABLE]

In other words, the extension $E(p)$ is a polynomial function with the same coefficients as $p$, but defined on the larger domain $K^{e}\cup K^{e}_{ext}$. We can finally address the issue of transferring the boundary conditions to the computational boundary. Let $\boldsymbol{q}_{h}$ be a polynomial approximation to $\boldsymbol{q}$. Then, for $\boldsymbol{x}\in\Gamma^{h}\cap e$ and $\overline{\boldsymbol{x}}\in\Gamma\cap K^{e}_{ext}$, the quantity

[TABLE]

will be used as an approximation to the boundary value $\psi(\boldsymbol{x})$. Note that the same formula can be used to define the approximation of the pointwise value of $\psi$ for any point $x\in K^{e}_{ext}$.

2.3 The discrete system

We can now present the form of (5) that will be discretized. Given a computational domain $\Omega^{h}$ with boundary $\Gamma^{h}$ and a regular embedded triangulation $\mathsf{T}^{h}$ as defined in Section 2.1, we look for functions $\psi$ and $\boldsymbol{q}$ defined on $\mathsf{T}^{h}$ and $\widehat{\psi}$ defined on $\mathcal{E}^{h}$ satisfying the system

[TABLE]

where $\varphi_{h}$ is given by (7) and the jump of the flux across two elements $K^{+},K^{-}$ with exterior normal vectors $\boldsymbol{n}^{+},\boldsymbol{n}^{-}$ along a shared edge $e$ is defined in standard fashion as

[TABLE]

. Note that for the remainder of this article, we will also use this notation for the jump of any scalar quantity $a$ across two elements, namely $[\![a]\!]:=a^{+}-a^{-}$. In (8c) and in what follows, quantities with the superscript $\;\widehat{\phantom{u}}\;$ must be understood as defined only on the skeleton of the mesh. The system (8) is the restatement of (5) as a collection of the local problems (8a) and (8b) satisfying the boundary conditions (8c) or (8e) on each element and “glued” together by the continuity condition on the flux (8d).

The introduction of the hybrid unknown $\widehat{\psi}$ as a global quantity encoding the boundary conditions allows to fully decouple the local problems if an appropriate numerical flux is chosen (more on this below). Once the hybrid variable has been determined, the local problems for $\psi$ and $\boldsymbol{q}$ can be solved independently. This process is akin to the well known static condensation technique used to decouple degrees of freedom on the edges/faces of an element from those on the interior which was first devised for Finite Elements [21] and mixed formulations through hybridization [22]. The connections between HDG and static condensation have been thoroughly discussed in [23].

When equation (8b) is posed weakly and discretized, one is faced with the choice of a numerical approximation for the normal flux across the element edges. Due to the fact that it allows to express the weak forms of (8a) and (8b) entirely in terms of local quantities and also that it allows great freedom of choice for the approximation spaces, the choice

[TABLE]

has become standard [23] and we will follow it in our discretization. Moreover, it is known that if the stabilization parameter $\tau>0$ remains of order $\mathcal{O}(1)$, the method achieves optimal convergence order. Therefore, for all the computations we will set it to be $\tau=1$.

2.4 The HDG discretization

The HDG method [24] yields piecewise polynomial approximations $(\boldsymbol{q}_{h},\psi_{h},\widehat{\psi}_{h})$ to the solutions of the weak formulation obtained by testing (8) with functions in the finite dimensional spaces

[TABLE]

where the space of polynomials of degree $k$ defined on the triangle $K$ is denoted by $\mathbb{P}_{k}(K)$, the product space of two copies of itself is given by $\mathbf{P}(K):=[\mathbb{P}_{k}(K)]^{2}$ and $\mathbb{P}_{k}(e)$ is the space of polynomials of degree $k$ defined on a given edge $e$. The $L^{2}$ inner products in these spaces are given by

[TABLE]

where, as is customary, $(\cdot,\cdot)_{K}$ and $\langle\cdot,\cdot\rangle_{\partial K}$ are the $L^{2}$ inner products on a single element $K$ and on its boundary $\partial K$ respectively.

Once the choice of trace of the numerical flux given by (9) has been introduced and the system has been tested with functions $(\boldsymbol{v},w,\mu)\in\boldsymbol{V}_{h}\times W_{h}\times M_{h}$, the weak form of the first order system (8) can be understood as consisting of two parts. The local equations

[TABLE]

that are satisfied by the hybrid unknown $\widehat{\psi}_{h}$ at the interior edges of the triangulation $\mathcal{E}^{\circ}_{h}$ (13c) and at the set of edges belonging to the boundary of the computational domain, $\mathcal{E}^{\partial}_{h}$ (13d). For each extended element $K_{ext}$ we can let $\overline{\boldsymbol{x}}\in\partial K_{ext}\cap\Gamma$ be the starting point of a transfer path and use the equation (7) together with the fact that $\psi$ satisfies homogeneous Dirichlet boundary conditions to express the transferred boundary value $\varphi_{h}$ appearing on (13d) as

[TABLE]

In the previous expression, $E(\boldsymbol{q}_{h})$ is the extrapolation of the polynomial $\boldsymbol{q}_{h}$ defined on $K^{e}$ to the neighboring exterior element $K_{ext}^{e}$ obtained by extending the domain of definition of $p$ from $K^{e}$ to $K^{e}_{ext}$ while keeping the same polynomial form.

2.5 The solution method

In this section we will first consider the source term $F$ to be independent of $\psi$; once the solution method has been described for this simple, linear case, we will come back to address the nonlinear case, in which $F=F(r,\psi)$. We will denote the basis functions of the approximation spaces $\boldsymbol{V}_{h}$, $W_{h}$, and $M_{h}$ respectively by $\boldsymbol{\phi}_{i},\phi_{i}$, and $\mu_{i}$, and define the Finite Element-style mass and convection matrices

[TABLE]

Moreover, if $\boldsymbol{q}_{h},\psi_{h},\widehat{\psi}_{h},\varphi_{h}$, and $F/r$ denote respectively the coefficient vectors of the system unknowns, the transferred boundary conditions, and the source term, we can write the two parts of the HDG system succinctly in matrix form as

[TABLE]

where the operators $\mathbf{I}_{e}^{\circ}$ and $\mathbf{I}_{e}^{\Gamma}$ above are the discrete counterparts of the restriction to the interior edges $\widehat{\psi}|_{\partial\text{T}^{h}\setminus\Gamma^{h}}$ and boundary edges $\widehat{\psi}|_{\Gamma^{h}}$, respectively. Now, following equation (14) the first term on the left hand side of (15t) comes from integrating the extrapolation $r\boldsymbol{E}(q_{h})$ along the transfer paths; it can therefore be expressed in the form

[TABLE]

where $\mathbf{Q}=\widetilde{\mathbf{Q}}\;\mathbf{I}_{K}^{\Gamma}$ is the composition of two discrete operators: 1) a restriction to the elements with at least one edge on the computational boundary, denoted as $\mathbf{I}_{K}^{\Gamma}$ , and 2) the combination of a line integral along the transfer paths and the inner product with the basis of $M_{h}$ defined on the skeleton of the mesh, which can be represented as

[TABLE]

At the implementation level, both integrals (the line integral and the inner product) involved in $\widetilde{\mathbf{Q}}$ are approximated by quadrature rules with matching orders of accuracy. Therefore the global system can be written in matrix form as

[TABLE]

where the top row is satisfied by the degrees of freedom of $\widehat{\psi}_{h}$ lying on the internal edges of the skeleton and the bottom row is satisfied by those on the edges corresponding to the boundary of the computational domain. Solving formally the linear system (15j) for $(\boldsymbol{q}_{h},\psi_{h})$ we obtain

[TABLE]

This expression for $(\boldsymbol{q}_{h},\psi_{h})$ can then be substituted into (16) and, from the resulting system, it follows that

[TABLE]

From this equation one can obtain $\widehat{\psi}_{h}$ and back-substitute in (17) to obtain $(\boldsymbol{q}_{h},\psi_{h})$. The relevance of the last two equations stems from the following observations. First, the matrices and the corresponding linear solves appearing on the right hand side of equation (17) are entirely in terms of local quantities, and can therefore be processed fully in parallel. Second, despite the fact that it involves a global unknown, the system (18) is sparse, for it includes only degrees of freedom associated to either the skeleton of the mesh or elements with at least one edge on the computational boundary. Moreover, the linear solves appearing on each side of equation (18) are the same as those in (17) and therefore have to be computed only once. The solution process can then be split into three steps:

Locally (i.e. in parallel) solve the systems

[TABLE]

appearing in equation (17) and store them. 2. 2.

Using the local vectors obtained in the first step, assemble the matrices on both sides of equation (18) and solve the resulting global system, thus recovering the hybrid unknown $\widehat{\psi}_{h}$. 3. 3.

Distribute the relevant parts of $\widehat{\psi}_{h}$ over local elements and use the local solvers obtained on the first step to recover $(\boldsymbol{q}_{h},\psi_{h})$ fully in parallel.

Accelerated Picard iterations

In order to deal with the non-linear nature of the Grad-Shafranov equation, we will resort to a simple, yet effective, iterative strategy consisting of accelerated Picard iterations. The standard Picard or fixed point iteration goes as follows. Given a guess for the solution $\psi^{n}$, the source term can be evaluated yielding a source $F:=F(r,\psi^{n})$ that is independent of the solution. The resulting system (15) is a linear problem that can be solved as described above yielding an update $\psi^{n+1}$. The source term is then updated by evaluation at the newly computed solution and the process is repeated iteratively until the relative difference between successive updates falls below a certain predetermined tolerance.

This strategy is simple to implement but may require a large number of iterations to converge for small values of the tolerance. However, the convergence rate for Picard iterations can be improved by means of a device known as Anderson acceleration [25]. Anderson’s idea is to improve convergence through the use of information from more than one previous iterate. This is achieved by defining the update $\psi_{h}^{n+1}$ to be an optimized convex linear combination of the solutions to (15) obtained on a predetermined number of previous iterations. The coefficients of the convex linear combination are chosen so that the difference between the solutions and the updates is minimized. This requires the storage of $m$ previous updates and $m$ previous solutions and the solution of a small $(m+1)\times(m+1)$ system at every iteration in order to determine the optimal coefficients.

Below we describe algorithmically the simplest form of the acceleration—which is the version implemented in our solver—but we refer the reader to the works by Kelly and Toth [26], and Walker and Ni [27], where the method is studied in detail. If we denote by $\epsilon$ a prescribed tolerance, by $\psi_{0}$ the initial input, by $\left(-\Delta^{*}\right)^{-1}$ the solution operator to (15) described above, and by $\psi_{h}$ the final, approximate solution to the non-linear problem then, in its simplest form, the acceleration algorithm that uses $m$ previous iterates can be described as follows:

2.6 Non-linear local post-processing

Following the idea introduced by Stenberg [28], once the approximations $\psi_{h}$ and $\boldsymbol{q}_{h}$ have been determined from the solution of (15), it is possible to define a locally post-processed function $\psi^{*}_{h}$ with enhanced convergence properties. There a several different ways of defining the post-processing, but in order to take advantage of the increased accuracy of the post processing as part of a residual estimator we will define $\psi_{h}^{*}$ to be the piecewise polynomial function satisfying

[TABLE]

Note that when $F$ is independent of $\psi$, this reduces to the case analyzed in [29], where it was shown that the solution to this auxiliary problem converges towards $\psi$ with order $k+2$ when $k\geq 1$. Numerical evidence suggests that the simpler post processing that arises if the terms involving $F$ in equation (19a) are dropped is also effective; even in the semi-linear case. However for our convergence analysis in [19] the effect of the non-linear source term needs to be considered and this leads to the non linear post processing above. The solution to this auxiliary problem will be used in the error estimator described in the next section.

3 The adaptive algorithm

In many situations of physical interest, the solution $\psi$ and its derivatives may vary rapidly in localized regions in $\Omega$ [16]. In such cases, adaptive mesh refinement is an effective way to minimize the number of degrees of freedom for a given target accuracy. Our adaptive strategy follows the standard “solve $\rightarrow$ estimate $\rightarrow$ mark $\rightarrow$ refine” iterative paradigm. Specifically, starting from an initial triangulation $\mathsf{T}^{h,0}$, the problem is solved and a suitable error estimator is computed using the obtained approximate solution. Based on the error estimator, a subset of the triangulation is marked for refinement. This generates a new triangulation $\mathsf{T}^{h,1}$ where the process can be started over until a predetermined number of cycles is reached or the estimator falls below a given threshold.

Our algorithm is based on theoretical work done by Cockburn and Zhang [30, 31], and Cockburn, Nochetto and Zhang [32]. In the first references, the authors proposed and studied the performance of a residual-based error estimator, and in the third one they were able to prove the convergence of the adaptive HDG method assuming that Dörfler’s marking criterion is used (we will return to this later). The problems studied in those cases were linear and the equations were posed in polygonal domains discretized with fitted triangulations. Our focus on the Grad-Shafranov equation poses additional challenges, namely the semi-linearity of the problem and the non-fitting nature of the computational domain—which in turn imposes the requirement that the distance between the boundaries $\Gamma$ and $\Gamma^{h}$ remains locally $\mathcal{O}(h_{loc}(\boldsymbol{y}))$. The non-linearity is dealt with through the accelerated iterative process described above. The refinement of the unfitted grid will require some additional care, as we discuss below.

3.1 A residual-type estimator

For an edge $e$ with length $h_{e}$ and a function $u$ defined on $e$ (or on a superset containing $e$) we will denote by $\|u\|_{e}$ its $L^{2}$ norm on $e$ (or the norm of its restriction to $e$). Considering $K$ to be a generic element of the triangulation $\mathsf{T}^{h}$, we will adopt the following local error estimator

[TABLE]

where $\psi^{*}_{h}$ is the post-processed numerical solution obtained by solving the local auxiliary problem (19). This estimator is based on the one proposed and analyzed by Cockburn and Zhang [30, 31] for linear elliptic equations posed in polygonal domains. The global error estimate is obtained by adding all the local contributions over the computational domain and is therefore defined as

[TABLE]

A detailed analysis of the estimator and its properties for semilinear problems like ours as well as possible improvements for it are the subject of a separate communication [19], but some intuitive understanding can be gained by expressing the estimator in the form $\eta^{2}=\eta_{1}^{2}+\eta_{2}^{2}+\eta_{3}^{2}+\eta_{4}^{2}+\eta_{5}^{2}$ , where

[TABLE]

and studying each term separately. The term $\eta_{1}$ corresponds to the local residual of the strong equation for the flux (8b). Similarly, it would be desirable to consider the residual $q_{h}-\frac{1}{r}\nabla\psi_{h}$ of the strong equation (8a) as part of the estimator. However, it would converge with reduced order due to the differentiation of the approximate solution $\psi_{h}$. To overcome this and achieve the desired order of convergence, the term $\eta_{2}$ makes use of the post processed solution $\psi^{*}_{h}$ instead, thus preserving the desired convergence order $k+1$. In this sense, the second term of the estimator is reminiscent of indicators based on gradient recovery, where an improved approximation of the gradient is obtained through post processing and is then used to estimate the error.

Finally, the edge terms $\eta_{3}$ and $\eta_{4}$ provide, respectively, a measure of the local loss of conformity of the solution and of its flux, by considering their jumps across element edges in the normal direction. $\eta_{5}$ estimates the rate of convergence of the hybrid variable and post-processed solution—restricted to the element boundaries—as approximations to the local trace.

3.2 Marking strategies

Given an initial triangulation $\mathsf{T}^{h,0}$, the discrete problem is solved and post-processed yielding the approximations $(\psi_{h,0},\boldsymbol{q}_{h,0},\widehat{\psi}_{h,0},\psi^{*}_{h,0})$ which are then used to compute the local error estimator (20). One now must choose the elements that will be marked for refinement based on the local values of $\eta_{K}$. Different marking strategies have been tried in the literature. Here we consider Dörfler’s criterion [33, 34] and the so-called maximum criterion [35]. In both cases, one must first choose a value of the marking parameter $\gamma\in[0,1]$ and then the elements

Either belonging to a minimal set $\mathcal{M}$ such that

[TABLE] 2. 2.

or for which the local estimate $\eta_{K}$ is such that

[TABLE]

are marked and subsequently refined. The choice of the value of the marking parameter $\gamma$ depends on the needs and constraints of the user. It is usually picked based on considerations such as memory availability, desired speed of the computation, etc. In the case of Dörfler’s marking, values of $\gamma$ closer to zero result in fewer elements being marked on each cycle, and very localized refinement; larger values of $\gamma$ tend to produce refinements that are more uniform. This qualitative behavior is reversed when maximum marking is used: refinement is localized for large values of $gamma$ and becomes uniform as the parameter approaches zero.

The convergence of the adaptive refinement loop on fixed polygonal domains was established by Cockburn, Nochetto and Zhang within the context of hybridizable discontinuous Galerkin methods for linear problems [32]; and assuming Dörfler’s method for marking. Regarding the maximum marking criterion, Morín, Siebert and Veeser [36] proved the convergence of the method for a wide class of linear problems discretized with Finite Elements. The analysis of this strategy applied to HDG is still an ongoing task, even for the linear case, but the strategy seems to be robust, as suggested by our numerical experiments.

3.3 Local mesh refinement

Once some elements of the triangulation have been marked by one of the criteria presented above, triangle refinement can be carried out through standard methods such as newest vertex bisection (NVB) [37, 38] or a red-green procedure [39]. In the standard setting where the computational domain $\Omega^{h}$ coincides with the domain of definition of the PDE, each of the refined meshes $\{\mathsf{T}^{h,n}\}_{n\geq 0}$ produced in such fashion will remain a triangulation of the original domain $\Omega_{h}$.

In the present situation however, the computational domain is in fact a strict subdomain of $\Omega$. Thus, if $\mathsf{T}^{h,0}$ is a triangulation of a fixed computational domain $\Omega_{h}$ built as described in Section 2.1, with every subsequent refinement step the mesh will drift farther away from satisfying the condition that the distance between the computational boundary $\partial\Omega_{h}$ and the actual boundary $\partial\Omega$ remains locally of the order of the triangle diameter $h_{loc}$, as depicted in Figure 3. As a result, the transfer procedure would not yield satisfactory results.

In order to avoid such situations, we propose a strategy to update the computational domain consistent with the mesh refinement in such a way that the local distance condition is always satisfied. The method, showed schematically in Figure 4, can be described as follows:

1. Starting from the computational domain $\Omega_{h,0}$, a pair of computational and companion meshes $\mathsf{T}^{h,0}$ and $\mathsf{T}^{h,0}_{c}$ are built following the process detailed in Section 2.1.

2. The problem is solved and the error is estimated on the mesh $\mathsf{T}^{h,0}$, which results in a list $\mathcal{M}$ of elements marked for refinement (green triangles in the second column of Figure 4).

3. The elements in the companion mesh $\mathsf{T}^{h,0}_{c}$ which correspond to those marked on the computational mesh in the previous step are marked for refinement. In addition, all the elements in the companion mesh which intersect the true boundary $\partial\Omega$ and share an edge with any triangle in $\mathcal{M}$ are marked for refinement as well (yellow triangles in the second column of Figure 4). This results in an augmented list $\mathcal{M}_{c}$ of elements in the companion mesh.

4. The companion mesh $\mathsf{T}^{h,0}$ is updated by performing triangle refinement on all elements in $\mathcal{M}_{c}$, yielding a temporary background mesh $\widetilde{\mathsf{T}}$ that completely contains $\Omega$, as depicted in the central column of Figure 4. Note that since the new elements may have smaller diameter, some of them may now in fact be completely contained in $\Omega$ even if the parent triangle was not. In a similar fashion, some of them may neither be contained in $\Omega$ nor intersect $\partial\Omega$.

5. A new computational domain $\Omega_{h,1}$ and its corresponding triangulation $\mathsf{T}^{h,1}$ are defined by selecting the triangles in $\widetilde{\mathsf{T}}$ that are completely contained in $\Omega$. Analogously, a new companion mesh $\mathsf{T}^{h,1}_{c}$ is defined by selecting all the elements of $\widetilde{\mathsf{T}}$ that are either completely contained in $\Omega$ or intersect the boundary $\partial\Omega$. The remaining triangles are discarded. The resulting level of refinement will then use the new computational and companion meshes and the process will continue until the predetermined stopping criterion is met.

Note that no computations are ever carried out using the companion grid, which is needed only to update the background triangulation in a way that gives rise to a refined computational domain $\Omega_{c,n}$ satisfying the separation condition $d(\boldsymbol{y},\partial\Omega)=\mathcal{O}(h_{loc}(\boldsymbol{y}))$ for every $\boldsymbol{y}\in\partial\Omega_{c,n}$. Moreover, since at every level the computational domain $\Omega_{h,n}$ is defined as a subdomain of an increasingly finer background triangulation, the algorithm yields a sequence of computational domains that effectively “exhaust” $\Omega$ as $n\to\infty$ so long as the local estimator $\eta_{K}$ remains nonzero on the elements having an edge on $\partial\Omega_{h,n}$. This is illustrated in Figure 5.

4 Numerical Experiments

We present five examples to showcase the performance of the adaptive algorithm. To establish the correct behavior of the error estimator and the refinement strategy we begin with a linear test case where the analytic solution is known. The chosen analytic equilibrium has the advantage of providing the desired geometric and parametric flexibility, but corresponds to solutions that vary smoothly across the confinement region and thus tend to favour a uniform refinement strategy. Nevertheless, the availability of an exact solution allows us to verify the overall behavior of the estimator.

More challenging and physically relevant situations arise when the source term of the equation is non-linear and has large gradients. The solutions for these cases develop features that would be hard to resolve accurately while simultaneously keeping a small number of elements and using the boundary transfer technique. It was precisely to address cases like these that we included adaptive refinement capabilities in our solver, and the last four examples in this section belong to this category. In all the experiments, the maximum criterion was used for marking.

4.1 A Solov’ev equilibrium

In order to test the performance of the error estimator and the adaptive algorithm, we start with a simple linear setting for which an exact solution is available. The example corresponds to a so-called Solov’ev profile [40] where the source term of the equation is taken to be of the form [41]

[TABLE]

The free parameter $A$ determines the ratio of plasma pressure to magnetic pressure in the equilibrium of interest. With such a source term, the Grad-Shafranov equation is linear and exact solutions can be constructed by imposing physical or geometrical constraints. The case that we will consider here corresponds to a geometry similar to that of the National Spherical Toroidal Experiment (NSTX) in the high beta regime.

The exact solution is of the form

[TABLE]

Graphs of the solution and its partial derivatives with these parameter values on the target geometry can be seen in Figure 6.

To test the behavior of the error estimator with respect to the true error, the equation was solved with polynomial basis with degrees $k=1,\ldots,4$ and uniform mesh refinement (i.e. for any subsequent levels $h_{n+1}=h_{n}/2$). The initial mesh diameter was $h=1.03$. The convergence plots for this experiment are shown in the top row of Figure 7, where it can be seen that the estimator accurately captures the qualitative behavior of the error. Moreover, as can be seen in Table 1, the rates of convergence of the numerical solutions $\psi_{h}$ and $\boldsymbol{q}_{h}$, as well as those of the error estimator $\eta$ and all its component terms $\eta_{1},\ldots,\eta_{5}$ are nearly optimal.

For comparison, five levels of the adaptive algorithm were ran on the same problem for polynomial degrees from 1 to 4 on the same initial grid. Marking was done using the maximum criterion with parameter $\gamma=0.3$: elements whose estimator is at least 30% of the maximum local estimate are marked. The convergence history can be seen in the bottom row of Figure 7. As can be seen in the same figure, for $k=4$ the number of degrees of freedom after four levels of uniform refinement is about the same order of magnitude as that of five levels of adaptive refinement, but the adaptive algorithm places most of the computational effort on the left side of the domain. Comparing with the plot of $\partial_{rr}\psi$ in Figure 6 it is clear that the refinement is focusing on the region where the magnitude of the second derivative in the horizontal direction peaks.

4.2 A pressure pedestal

The following example features a pressure profile that remains almost flat throughout the confinement region and drops abruptly in the vicinity of the boundary:

[TABLE]

with $c_{1}=0.8$, $c_{2}=0.2$, and $\sigma^{2}=0.005$. This kind of profile is quite frequent in magnetic confinement fusion experiments, where the narrow region of fast decrease of the pressure is known as a pressure pedestal. If the equilibrium is assumed to be neither paramagnetic nor diamagnetic, $g(\psi)=constant$, and the source term of the Grad-Shafranov equation is

[TABLE]

which has very strong gradients close to the edge of the confinement region since $\sigma$ is small. Figure 8 shows both the pressure and the source profiles in an ITER-like geometry with an magnetic X-point (as described in [41]). As can be seen from the cross sections at constant values of $z$ in the same figure, the source has large gradients close to the boundaries, especially in the “outer” region.

The equation was solved using (24) as the source term and the latter parameter values. Figure 9 shows the post-processed numerical solution obtained when the polynomial basis was chosen to have degree four, and six levels of refinement were used with marking parameter $\gamma=0.3$. The computational mesh had initial mesh parameter $h=1.71\times 10^{-1}$, whereas the final mesh, displayed in the first block of the figure, consists of 636 elements with maximum diameter $h_{max}=1.71\times 10^{-1}$ and minimum diameter $h_{min}=1.07\times 10^{-2}$. Our adaptive refinement scheme clearly focuses degrees of freedom in the region corresponding to the pressure pedestal, where $\psi$ and its derivatives vary strongly.

4.3 A transport barrier

In magnetic confinement fusion experiments, large pressure gradients may also be observed closer to the core of the discharge, and correspond to internal transport barriers [42]. To model such situations, we consider a pressure profile of the form

[TABLE]

where erf$(\cdot)$ is the error function, $a$ and $b$ are natural numbers, and the parameters $H$ and $s$ control the height and the steepness of the barrier respectively. The parameter $\psi_{0}$ gives the location of the transport barrier with respect to a common normalization where $\psi\in[0,1]$ [43, 17]. In Figure 10 we present the behavior of the pressure as a function of $\psi$ for different parameter values, and thus verify that we are able to model experimentally relevant situations [42, Figure 3].

Figure 11 shows the numerical solution for a transport barrier with parameter values $H=0.5$, $s=40$, $\psi_{0}=0.3$, $a=4$, and $b=2$, corresponding to the steepest barrier depicted in Figure 10. Plots of the pressure distribution and the source term in the confinement geometry are shown side by side. The geometry for this experiment corresponds to an up-down symmetric ITER-like configuration with two magnetic X-points [41]. The solution displayed in the figure was computed adaptively using piecewise cubic polynomials. The initial mesh consisted of 129 elements, and the final grid consisted of 514 elements with a ratio $h_{max}=8.97\times 10^{-2}$ and $h_{min}=7.9\times 10^{-3}$. Both meshes can be seen in Figure 12.

Figure 12 shows the convergence of the global error estimator for basis functions of degrees varying from 1 to 4. All the experiments started on the initial mesh shown in Figure 12, which subsequently underwent six levels of refinement with maximum marking using $\gamma=0.5$ (i.e. only those elements whose local contribution is at least 50% of the maximum value of $\eta_{k}$ are refined at every level). Note that at the finest level of refinement all the examples in the figure have roughly the same number of degrees of freedom; however, as is expected, higher order discretizations result in coarser grids (i.e. consisting of fewer elements) and considerably smaller error.

4.4 A current hole

Another challenging configuration of physical relevance is the family of equilibria with a so-called “current hole”. This name refers to equilibria for which there is an extended region in the plasma core where the toroidal current $J_{\phi}$ is nearly zero [44, 45, 46, 47, 48]. The “near absence” of the current (thus the name “hole”) corresponds to a core region with almost constant pressure. Such equilibria have been observed experimentally in several tokamaks, and are found to be remarkably robust [48]. They are likely to naturally occur in future large scale tokamaks in fully non-inductive current drive operation. They have therefore gathered significant interest in recent years [48] (and references therein).

The following source term, adapted from that of the pressure pedestal, can give rise to such an equilibrium

[TABLE]

Choosing the values $c_{1}=0.4,c_{2}=0.1,c3=-18,c_{4}=10\pi$ as well as $\sigma^{2}_{1}=5\times 10^{-3}$ and $\sigma^{2}_{2}=3\times 10^{-3}$, the spatial distribution of the source term will be as shown in Figure 13 for an up-down symmetric D-shaped geometry with ITER-like parameters (as described in [16]). Recalling the relationship $\mu_{0}J_{\phi}=F(r,\psi)/r$, we see that the current drops to close to zero in the core of the confinement region and has sharp peaks near the boundary. As can be seen in the figure, there is a sharp contrast between the almost constant behavior of both the source and the solution in the central region and the large gradients at the edge. Consistently, the estimator focuses the computational effort on the edges and keeps a relatively coarse mesh in the core. The computation was ran over six levels of refinement with marking parameter $\gamma=0.3$; for a polynomial basis of degree four the resulting final mesh shown in the figure consisted of 3310 elements where $h_{min}=5.1\times 10^{-3}$ and $h_{max}=8.16\times 10^{-2}$.

Figure 14 shows the partial derivatives of $\psi_{h}$ obtained directly from the components of the flux $\boldsymbol{q}_{h}$ and their cross sections for different values of $z$. The second derivatives were computed by differentiating the local polynomial approximants of $\partial_{r}\psi_{h}$ and $\partial_{z}\psi_{h}$, which introduces additional error. However, the high polynomial order in combination with the focused mesh refinement imply that, even with the expected deterioration, the approximation remains within an acceptable range.

4.5 An internal layer

Perhaps one of the most desirable features in an adaptive scheme is the ability to automatically detect localized features in the solution and to refine the computation locally in order to resolve them accurately. The pressure pedestal can be turned into a more challenging benchmark along these lines if the source term is modified to

[TABLE]

This source term is not physically relevant for magnetic confinement fusion applications, because it cannot be cast in the canonical form of the source in (1) due to the explicit appearance of the coordinate $r$ in the argument of the last exponential. Nevertheless, it represents a good benchmarking problem to test for the detection of internal layers. As can be seen in Figure 15 (left and center left) the source presents an internal layer that changes abruptly in addition to the large gradients at both edges of the confinement region. In the figure and the numerical experiment the constants were taken to be $c_{1}=0.8,c_{2}=0.2,c_{3}=15,\sigma^{2}_{1}=5\times 10^{-3},$ and $\sigma^{2}_{2}=7.5\times 10^{-4}$.

The simulation parameters were as in the example with the pressure pedestal: the same ITER-like geometry with an x-point, the same starting grid, polynomial basis of degrees one to four and six levels of refinement with $\gamma=0.3$. As can be seen in the central panel of Figure 15, the estimator successfully detects the development of internal features in the solution, and concentrates the refinement in that region. The final grid consists of only 601 elements with maximum mesh diameter $h_{max}=1.21\times 10^{-1}$ and minimum mesh diameter $h_{min}=2.14\times 10^{-2}$.

The post-processed numerical solution and cross sections at different heights are depicted on the right end of Figure 15. The sharp change in the slope of the solution drives the interior refinement thus enabling the accurate approximation of the step-like behavior of the derivative in the horizontal direction (Figure 16). The approximate first derivatives are shown in Figure 16 along with a plot of the convergence history of the global error estimator.

5 Conclusion

The solver for fixed-boundary Grad-Shafranov equilibria based on the hybridizable discontinuous Galerkin method we presented has several attractive features beyond the high order approximation properties for $\psi$ and its partial derivatives. The use of an HDG framework provides the code with a robust and reliable method that is naturally suited for parallelization and addresses the issue of unnecessarily large numbers of degrees of freedom—usually associated with discontinuous Galerkin discretizations—through hybridization.

The use of a polygonal subdomain to carry out the computations combined with the transfer algorithm to impose boundary conditions “at a distance” allows for a simple, yet highly accurate, treatment of curved boundaries without having to resort to more complicated techniques like isoparametric mappings. It also enables a unified treatment of both smooth geometries and those with corners, corresponding to magnetic X-points. Moreover, in applications where the geometry of the confinement region needs to be updated, this technique provides the additional benefit of avoiding the need for constant updating of a fitted mesh.

The addition of a local error estimator and adaptive mesh refinement allows for focused computational efforts. As the numerical experiments show, this feature combined with the updating of the computational domain and the improved geometric approximation of the physical domain as the refinement progresses allows for efficient detection of relevant physical effects near the edge (as is the case for equilibria with a pressure pedestal) or the resolution of highly localized internal structures (as is the case for equilibia with an internal transport barrier).

6 Acknowledgements

The computational implementation of the algorithm described in this paper benefited greatly from the detailed explanations and code templates for HDG and adaptive refinement provided respectively by Fu, Gatica and Sayas [49], and Funken, Praetorius and Wissgott [50]. Finally, sampling of the confinement regions from the analytic expressions given in [41] was done using chebfun [51].

The authors are deeply grateful to Wrick Sengupta (NYU) and Georg Stadler (NYU) for their valuable insights on the physical and mathematical aspects of the problem. They also thank François Waelbroeck (UT-Austin) for suggesting the current hole problem as a benchmarking test. Antoine J. Cerfon. and Tonatiuh Sánchez-Vizuet were partially funded by the US Department of Energy. Grant No. DE-FG02-86ER53233. Manuel E. Solano was partially funded by CONICYT–Chile through FONDECYT project No. 1160320 and by Project AFB170001 of the PIA Program: Concurso Apoyo a Centros Cientificos y Tecnologicos de Excelencia con Financiamiento Basal.