An adaptive discontinuous Petrov-Galerkin method for the Grad-Shafranov equation
AAN ADAPTIVE DISCONTINUOUS PETROV-GALERKIN METHODFOR THE GRAD-SHAFRANOV EQUATION ∗ ZHICHAO PENG † , QI TANG ‡ , AND
XIAN-ZHU TANG § Abstract.
In this work, we propose and develop an arbitrary-order adaptive discontinuousPetrov-Galerkin (DPG) method for the nonlinear Grad-Shafranov equation. An ultraweak formula-tion of the DPG scheme for the equation is given based on a minimal residual method. The DPGscheme has the advantage of providing more accurate gradients compared to conventional finiteelement methods, which is desired for numerical solutions to the Grad-Shafranov equation. The nu-merical scheme is augmented with an adaptive mesh refinement approach, and a criterion based onthe residual norm in the minimal residual method is developed to achieve dynamic refinement. Non-linear solvers for the resulting system are explored and a Picard iteration with Anderson accelerationis found to be efficient to solve the system. Finally, the proposed algorithm is implemented in parallelon MFEM using a domain-decomposition approach, and our implementation is general, supportingarbitrary order of accuracy and general meshes. Numerical results are presented to demonstrate theefficiency and accuracy of the proposed algorithm.
Key words. discontinuous Petrov-Galerkin method, adaptive mesh refinement, high-order,Grad-Shafranov
AMS subject classifications.
1. Introduction.
The magnetohydrodynamic (MHD) equilibrium is critical inmany applications of plasma systems, such as magnetic confinement fusion. Steady-state fusion reactor operation requires that a MHD equilibrium is reached and sus-tained. In computational plasma physics, an efficient and accurate solver for the MHDequilibrium serves as the basis for linear stablity analysis and nonlinear simulationsof plasma transport and off-normal events such as tokamak disruptions. For example,the confining magnetic fields computed from a MHD equilibrium solver can be usedas a starting point for the relativistic Fokker-Planck equation in the runaway electronstudy [23], or a starting point to evolve a nonlinear MHD simulation [27]. For moredetails on the importance to have accurate and efficient equilibrium solvers, see thediscussions in [32] for instance.The axisymmetric MHD equilibrium is governed by the Grad-Shafranov equation,which is an elliptic equation with a nonlinear source term [22, 50]. There are two typesof equilibrium problems, the so-called fixed boundary and free boundary equilibria. Inthe fixed boundary equilibrium problem, it is assumed that the separatrix of the MHDequilibrium is known and a Dirichlet boundary condition is imposed at the separatrixfor the Grad-Shafranov equation. (In practice, a closed flux surface just inside theseparatrix, the so-called q surface, is used for the computational boundary.) Thefree boundary problem considers the case of an unbounded domain when separatrixis not known and is in fact to be determined self-consistently from both the plasmacurrent and electrical currents in the external field coils. Thus, a free boundary solver ∗ Submitted to the editors DATE.
Funding:
This work was supported by the U.S. Department of Energy through the FusionTheory Program of the Office of Fusion Energy Sciences, and the Tokamak Disruption Simulation(TDS) SciDAC partnerships between the Office of Fusion Energy Science and the Office of AdvancedScientific Computing. † Department of Mathematical Sciences, Rensselaer Polytechnic Institute, Troy, NY 12180, USA.([email protected]). ‡ Los Alamos National Laboratory, Los Alamos, NM 87545, USA. ([email protected]). § Los Alamos National Laboratory, Los Alamos, NM 87545, USA. ([email protected]).1 a r X i v : . [ m a t h . NA ] J u l Z. PENG, Q. TANG, AND X.-Z. TANG needs extra capabilities, such as a boundary integral term that addresses the effect ofa far-field boundary and numerical routines to locate the separatrix in a given trialsolution. As an example, see [24] for a finite-element-based free boundary solver.In this work, we focus on the fixed boundary problem and develop a fast and high-order finite-element-based solver on general meshes. The fixed boundary solver is anecessary building block for a free boundary solver.In recent years, the fixed boundary solvers have drawn significant attentions due toits importance in magnetic confinement fusion. Many numerical schemes have beenexplored and developed. These include hybridizable discontinuous Galerkin meth-ods [48, 49], spectral elements [26, 38], boundary integral approaches [39, 32], Hermitefinite element [28, 33] and so on. Some of them, such as [49], also considered adap-tive mesh refinement (AMR) approaches to minimize numerical errors and improvethe efficiency of the algorithm. The derivative of the solutions determines physicallyimportant quantities such as magnetic fields for computing the charged particle tra-jectory. Hence, other than efficiency and accuracy of the algorithms, many of theprevious works also focus on minimizing errors in the derivatives of the solutions,including using bicubic elements or introducing auxiliary variables for its derivatives.Constructing accurate derivatives and effective AMR will be the two primary objec-tives of the current work.Among many available high-order schemes, we are interested in using discontin-uous Petrov Galerkin (DPG) methods for the Grad-Shafranov equation. The DPGmethod proposed in [13, 14] enjoys the following properties that are particularly at-tractive for our purpose. First, with an ultra-weak formulation and suitable choice ofthe discrete test space, the DPG method can provide high order accurate approxima-tion to both solutions and the first derivatives [21]. Second, the DPG method providesa natural built-in error estimator for adaptivity based on a numerical residual [17, 8],while a standard AMR method relies on the calculation of local numerical fluxes orphysical features like high order derivatives.Next we briefly review the DPG methodology. More details on the formulationscan be found in Section 3.1. The DPG method is a Petrov-Galerkin finite elementmethod, which uses on-the-fly computed optimal test functions [13, 14]. Discontin-uous or broken test space is used so that these optimal test functions can be locallycomputed on each element. Besides the perspective of optimal functions, the DPGmethod can be equivalently seen as a minimal residual method [18, 16, 30, 11, 44]. TheDPG method preserves a discrete version of the inf-sup condition and provides thebest error approximation for an energy norm defined in the trial space [18]. Anotherattractive feature of the method is that it provides a natural built-in error estima-tor [17, 8]. AMR and hp -adaptivity algorithms [17, 11, 44, 41] are designed based onthis estimator. DPG methods for nonlinear problems have also been considered. Forinstance, in [12, 11, 44], the authors first linearize nonlinear problems and then applythe DPG method to the linearized problems. In [35], a Picard iteration is applied tosolve DPG method for a nonlinear optical problem. In [34], from the perspective ofminimal residual method, a hybridized DPG method is designed and directly appliedto nonlinear fluid problems, and a Newton method is used for the nonlinear solver.In [6], a partial-differential-equation-constrained optimization method is consideredas a nonlinear solver for the DPG method. Recently, a posteriori error analysis is alsogeneralized to the nonlinear problem in [7].In this work, we would like to take full advantage of the DPG methodology asdescribed above and develop an efficient and accurate solver for the Grad-Shafranovequation, which produces not only high-order accurate solution but also high-order DAPTIVE DPG FOR GRAD-SHAFRANOV mfem.org ), and the solver supportsarbitrary order accuracy and general meshes for complex geometry.The remainder of the paper is structured as follows. The MHD equilibrium andGrad-Shafranov equation are briefly described in Section 2. The formulation of theDPG scheme is given in Section 3. Our discussion starts with a review of the abstractDPG formulation based on minimal residual method. The ultraweak formulation isthen described for the Grad-Shafranov equation. The details on discretizations suchas its matrix-vector form follow. In Section 4, we focus on efficient nonlinear solversfor the resulting matrix-vector form of the DPG scheme. In Section 5, we discussadaptive mesh refinement, followed by details of the implementation in Section 6.Finally, numerical results are presented in Section 7.
2. Governing equation.
In this section, we briefly introduce the MHD equi-librium and Grad-Shafranov equation, and then define the fixed-boundary problemthat is solved by the DPG scheme in this work.A MHD equilibrium for magnetically confined plasma means the force balancing j × B = ∇ p, where B is the magnetic field, j is the plasma current density that satisfies Ampere’slaw µ j = ∇ × B , with µ the magnetic permeability, and p the plasma pressure. If further assumingthe problem is axisymmetric (e.g., the equilibria in tokamaks) and considering theequilibrium in ( r, z )-coordinate, the Grad-Shafranov equation for the magnetic fluxfunction ψ can be written as r ∂∂r (cid:18) r ∂ψ∂r (cid:19) + ∂ ψ∂z = − µ r dpdψ − I dIdψ , where I is a function of ψ and it is associated with the toroidal part of the magneticfield. The magnetic field satisfies B = 1 r ∇ ψ × ˆ e θ + Ir ˆ e θ Note that the quantities of interest in practice are the derivatives of ψ (correspondingto the magnetic field) and its second derivatives (corresponding to the current). Formore details on the Grad-Shafranov equation and plasma equilibrium, the readers arereferred to plasma textbooks such as [29].Define a gradient operator (cid:101) ∇ = ( ∂ r , ∂ z ) T . In this work, we consider the Grad-Shafranov equation with the fixed boundary condition and the problem is rewrittenas (cid:101) ∇ · (cid:18) r (cid:101) ∇ ψ (cid:19) = − F ( r, z, ψ ) r , x ∈ Ω ,ψ = ψ D , x ∈ ∂ Ω , (2.1) Z. PENG, Q. TANG, AND X.-Z. TANG where Ω ⊂ R is the physical domain with a Lipschitz boundary ∂ Ω. In tokamaks,the physical domain Ω corresponds to the cross section of the device. The sourceterm F ( r, z, ψ ) depends on p and I , and in practice both of them are functions of ψ ,which are from either experimental measurements or theoretical design. Therefore,the problem is nonlinear through the source term and constructing efficient nonlinearsolvers is one of the main focuses of the current work.
3. Formulation of the DPG method.
We begin with a quick review of theformulation of the DPG method as a minimal residual method for a general abstractnonlinear problem. Based on this abstract formulation, we will show some details ofan ultraweak DPG formulation for the Grad-Shafranov equation.
Consider anabstract weak formulation for a general nonlinear problem. Let U denote the trialspace and V the test space, both of which are Hilbert spaces. The weak formulationfor a nonlinear problem is to find u ∈ U such that, b N ( u, v ) = l ( v ) , ∀ v ∈ V, (3.1)where l ( · ) : V → R is a linear form and b N ( u, v ) : U × V → R is nonlinear for u ∈ U while linear for v ∈ V . A residual operator can be defined as r ( u, v ) = b N ( u, v ) − l ( v ) . (3.2)Let V (cid:48) denote the dual space of V . Since b N ( u, v ) is linear for v , based on the Rieszrepresentation theory, there exists a linear operator B : U → V (cid:48) such that ∀ u ∈ U , < B u, v > V (cid:48) × V = b N ( u, v ), ∀ v ∈ V . It is also convenient to define the inner products( · , · ) U and ( · , · ) V as the inner products in U and V , respectively.As pointed out in [34, 6, 30, 7] for nonlinear problems, the DPG method is equiv-alent to a minimal residual method. In this work, we focus on this interpretationof the DPG method in both the discussion and implementation. Here, we brieflyreview the minimal residual interpretation, and more discussions can be found in[34, 18, 11, 43, 16, 30]. Suppose U h ⊂ U and V h ⊂ V to be the finite dimensionaldiscrete trial and test spaces. The DPG method can be defined as looking for u h ∈ U h such that u h = arg min w h ∈ U h ||B ω h − l || V (cid:48) = arg min w h ∈ U h (cid:32) sup v h ∈ V h , v h (cid:54) =0 | b N ( w h , v h ) − l ( v h ) | || v h || V (cid:33) = arg min w h ∈ U h (cid:32) sup v h ∈ V h , v h (cid:54) =0 | r ( w h , v h ) | || v h || V (cid:33) . (3.3)Let { e v i } Ni =1 be a basis of the discrete test space V h and e v := ( e v , . . . , e v N ) T , thenfor ∀ v h ∈ V , it can be expanded as v h = v T e v = (cid:80) Ni =1 v i e v i . Since both l ( v ) and b N ( u, v ) are linear with respect to the argument v , we have r ( w h , v h ) = r (cid:32) w h , N (cid:88) i =1 v i e v i (cid:33) = N (cid:88) i =1 v i r ( w h , e v i ) = v T r ( w h ) , DAPTIVE DPG FOR GRAD-SHAFRANOV r ( w h ) := ( r ( w h , e v ) , . . . , r ( w h , e v N )) T . Then, we note that | r ( w h , v h ) | || v h || V = | v T r ( w h ) | v T G v , (3.4)where the matrix G ∈ R N × N is defined as G ij := (cid:0) e v i , e v j (cid:1) V . (3.5) G is often called a Gram matrix, and here it can be interpreted as a mass matrixcorresponding to the inner product ( · , · ) V .The Gateaux derivative of (3.4) with respect to v is2( v T G v ) (cid:104) ( v T r ( w h ))( v T G v ) r ( w h ) − ( v T r ( w h )) G v (cid:105) . (3.6)It is easy to see when v satisfies v = v T G vv T r ( w h ) G − r ( w h ) , the Gateaux derivative is 0, and the supremum of | r ( w h ,v h ) | || v h || V is obtained assup v h ∈ V h , v h (cid:54) =0 | r ( w h , v h ) | || v h || V = r ( w h ) T G − r ( w h ) . (3.7)Therefore, the DPG method (3.3) is equivalent to find u h = arg min w h ∈ U h (cid:16) r ( w h ) T G − r ( w h ) (cid:17) . (3.8)Taking the Gateaux derivative of (3.8) gives J T ( u h ) G − r ( u h ) = 0 , (3.9)where J ( u h ) is the Jacobian matrix of the residual r ( u h ) with repect to u h . Finally,the DPG method based on the minimal residual method is given by (3.9), which willbe the focus of the current work. More details for the minimal residual interpretationof the DPG methodology can be found in [34, 30].The original DPG scheme was defined by computing optimal test functions onthe fly, see [14] for instance. The DPG scheme in [14] is equivalent to the form (3.9)derived from the minimal residual method. We refer the readers to [34] for details ofits equivalence proof. In summary, the methodology of the DPG method uses special“optimal” test functions, which preserves a discrete inf-sup condition with a discreteinf-sup constant of the same order as the original continuous inf-sup constant andguarantees the “optimal convergence” results, see [14, 18] for instance.Finally, consider the special case of linear problems. Let { e u i } Mi =1 be the basisof the discrete trial space U h and e u := ( e u , . . . , e u M ) T . For ∀ u h ∈ U h , it can beexpanded as u h = u T e u . For a linear problem when b N ( u, v ) becomes a bilinear form,one can verify that (3.9) becomes B T G − B u = B T G − (cid:96), where the matrix B ∈ R N × M defined by B ij = b ( e u j , e v i ) and (cid:96) is a constant vectordetermined by operator l ( · ). Hence, for linear problems, the DPG method alwaysresults in a symmetric positive definite linear system, while nonlinear problems donot have such a property, which will be discussed in Section 4. Z. PENG, Q. TANG, AND X.-Z. TANG
We present the details of the DPG formulation for the problem (2.1). The motivationto use the DPG scheme is to obtain a more accurate approximation to the gradientof ψ , and we therefore consider the ultraweak formulation to discretize (2.1). Definean auxiliary variable of the vector q := − (cid:101) ∇ ψr , and rewrite (2.1) into its first order form, r q + (cid:101) ∇ ψ = 0 , x ∈ Ω , (3.10a) (cid:101) ∇ · q = 1 r F ( r, z, ψ ) , x ∈ Ω , (3.10b) ψ = ψ D , x ∈ ∂ Ω . (3.10c)Let Ω h = { K } be a partition of the physical domain Ω, where K ∈ Ω h are disjointelements. Let ∂ Ω h = { ∂K, K ∈ T h } denote its skeleton (edge) and Γ h = ∂ Ω h ∩ ∂ Ω be the set of the edges on the physical boundary. Define ( · , · ) K and (cid:104)· , ·(cid:105) ∂K asthe standard L inner product of L ( K ) and L ( ∂K ), respectively. Let ( · , · ) Ω h = (cid:80) K ∈ Ω h ( · , · ) K , (cid:104)· , ·(cid:105) ∂ Ω h = (cid:80) ∂K ∈ ∂ Ω h (cid:104)· , ·(cid:105) ∂K and (cid:104)· , ·(cid:105) Γ h = (cid:80) ∂K ∈ Γ h (cid:104)· , ·(cid:105) ∂K .We generalize the DPG method for the linear Poisson problem in [21] to thenonlinear Grad-Shafranov equation. Define the trace spaces H − ( ∂ Ω h ) := { (cid:98) q n : ∃ q ∈ H (div; Ω) such that (cid:98) q n = q · n | ∂ K , ∀ K ∈ Ω h } ,H ( ∂ Ω h ) := { (cid:98) ψ : ∃ ψ ∈ H (Ω) such that (cid:98) ψ = ψ | ∂ K , ∀ K ∈ Ω h } , and the broken Sobolev spaces H (div , Ω h ) := { φ ∈ L (Ω h ) : φ | K ∈ H (div , K ) , ∀ K ∈ Ω h } ,H (Ω h ) := { τ ∈ L (Ω h ) : τ | K ∈ H ( K ) , ∀ K ∈ Ω h } . The trial space U and the test space V are chosen as U = (cid:0) L (Ω) (cid:1) × L (Ω) × H − ( ∂ Ω h ) × H ( ∂ Ω h ) ,V = H (div , Ω h ) × H (Ω h ) . The ultraweak form associated with (3.10) is to find u = ( q , ψ, (cid:98) q n , (cid:98) ψ ) ∈ U such thatfor ∀ v = ( φ , τ ) ∈ V ( r q , φ ) Ω h − (cid:16) ψ, (cid:101) ∇ · φ (cid:17) Ω h + (cid:104) (cid:98) ψ, n · φ (cid:105) ∂ Ω h = 0 , (3.11a) − ( q , (cid:101) ∇ τ ) Ω h + (cid:104) (cid:98) q n , τ (cid:105) Ω h = (cid:18) F ( r, z, ψ ) r , τ (cid:19) Ω h , (3.11b) (cid:104) (cid:98) ψ, τ (cid:105) Γ h = (cid:104) ψ D , τ (cid:105) Γ h . (3.11c)The source term F ( r, z, ψ ) can be rewritten as the summation of a nonlinear part anda linear part F ( r, z, ψ ) = F N ( r, z, ψ ) + F L ( r, z ). Based on the weak form (3.11), we DAPTIVE DPG FOR GRAD-SHAFRANOV l ( v ) = (cid:18) F L ( r, z ) r , τ (cid:19) Ω h , (3.12a) b N ( u , v ) = ( r q , φ ) Ω h − (cid:16) ψ, (cid:101) ∇ · φ (cid:17) Ω h + (cid:104) (cid:98) ψ, n · φ (cid:105) ∂ Ω h − ( q , (cid:101) ∇ τ ) Ω h + (cid:104) (cid:98) q n , τ (cid:105) Ω h − (cid:18) F N ( r, z, ψ ) r , τ (cid:19) Ω h . (3.12b)In addition, we define the test norm || · || V as || v || V := || ( φ , τ ) || V = || φ || + || (cid:101) ∇ · φ || + || τ || + || (cid:101) ∇ τ || , (3.13)and let ( · , · ) V be its corresponding inner product.Finally, we need to determine the discrete trial space U h and the discrete testspace V h . Let P k ( K ) and P k ( ∂K ) be the space of polynomials with degree at most k on the element K and its edge ∂K . The discrete trial space U h and the discrete testspace V h are respectively chosen as U kh = (cid:110) u h = ( q h , ψ h , (cid:98) q n,h , (cid:98) ψ h ) : q h | K ∈ (cid:0) P k ( K ) (cid:1) , ψ h | K ∈ P k ( K ) , (cid:98) q n,h | ∂K ∈ P k ( ∂K ) ∩ H − ( ∂ Ω h ) , (cid:98) ψ h | ∂K ∈ P k +1 ( ∂K ) ∩ H ( ∂ Ω h ) (cid:111) ,V k,sh = (cid:110) v h = ( φ h , τ h ) : φ h | K ∈ (cid:0) P k + s ( K ) (cid:1) , τ h | K ∈ P k + s ( K ) (cid:111) , s ≥ . According to [8], the key to prove the convergence of the DPG scheme is to definea Fortin-type operator Π : V → V h = V k,sh such that b N ( w h , ( I − Π) v ) = 0 for ∀ w h ∈ U h , v ∈ V . As proved in [21], s ≥ s ≥ U h , V h , (3.12) and (3.13),we are ready to determine G , J ( u h ) and r ( u h ) in (3.9) and define our DPG scheme. Remark || · || V is not unique. An attractivealternative is the adjoint graph norm [52, 18, 40]: || v || V, adjoint := || ( φ , τ ) || V, adjoint = || r φ − (cid:101) ∇ τ || + || (cid:101) ∇ · φ || + || φ || + || τ || , (3.14)which delivers a quasi-optimal error estimate in L sense. In [15, 21], our currentchoice ||·|| V is also proved to provide a quasi-optimal error estimate for linear ellipticalproblems. In this section, the details ofthe DPG scheme are presented as a matrix-vector form. The following notations willbe used throughout this section. Let U h ⊃ { e u j } Mj =1 := { e q j q } M q j q =1 × { e ψj ψ } M ψ j ψ =1 × { e (cid:99) q n j (cid:100) qn } M (cid:100) qn j (cid:100) qn =1 × { e (cid:98) ψj (cid:98) ψ } M (cid:98) ψ j (cid:98) ψ =1 and V h ⊃ { e v i } Ni =1 := { e φ i φ } N φ i φ =1 × { e τi τ } N τ i τ =1 Z. PENG, Q. TANG, AND X.-Z. TANG be the bases of the discrete trial space U h and the test space V h , respectively. Notethat e q j and e φ i are vector-valued, while other basis functions are scalar-valued. Define e ξ u := ( e ξ , . . . , e ξM ξ ) T and e η v := ( e η , . . . , e ηN η ) T , where ξ = q , ψ, (cid:98) q n , (cid:98) ψ and η = φ , η .Then, q h , ψ h , (cid:98) q n,h and (cid:98) ψ h can be rewritten as q h = Q T e qu , ψ h = Ψ T e ψ u , (cid:98) q n,h = (cid:99) Q n T e (cid:99) q n u , (cid:98) ψ h = (cid:98) Ψ T e (cid:98) ψ u . Define U := ( Q T , Ψ T , (cid:99) Q n T , (cid:98) Ψ T ) T , and on each element Q contains two componentscorresponding to the vector q . Now, we are ready to rewrite the ultraweak formula-tion (3.11) into its matrix-vector form and obtain the residual as r ( u h ) = b N ( u h , v h ) − l ( v h )= B L U − B N ( U ) − F L = (cid:18) M r (cid:101) ∇ h T (cid:98) ψ div h T (cid:99) q n (cid:19) QΨ (cid:99) Q n (cid:98) Ψ − (cid:18) N ( Ψ ) (cid:19) − (cid:18) L (cid:19) , (3.15)where the block matrices M r , (cid:101) ∇ h , T (cid:98) ψ , div h and T (cid:99) q n are defined as( M r ) ij := (cid:16) e q j , e φ i (cid:17) Ω h , ( (cid:101) ∇ h ) ij := − (cid:16) e ψj , (cid:101) ∇ · e φ i (cid:17) Ω h , ( T (cid:98) ψ ) ij := (cid:104) e (cid:98) ψj , n · e φ i (cid:105) ∂ Ω h , (div h ) ij := − (cid:16) e q j , (cid:101) ∇ e τi (cid:17) Ω h , ( T (cid:99) q n ) ij := (cid:104) e (cid:99) q n j , e τi (cid:105) ∂ Ω h , and the vectors N ( Ψ ) and L are defined as( N ( Ψ )) i := ( F N ( r, z, ψ h ) /r, e τi ) Ω h , L i := ( F L ( r, z ) /r, e τi ) Ω h . The Jacobian matrix of r ( u h ) can be therefore derived as J ( u h ) = (cid:18) M r (cid:101) ∇ h T (cid:98) ψ div h − D N ( ψ h ) T (cid:99) q n (cid:19) (3.16)where the block matrix D N ( ψ h ) is defined as ( D N ( ψ h )) ij := (cid:16) ∂ ( F N ( r,z,ψ h ) /r ) ∂ψ e ψj , e τi (cid:17) Ω h and ∂ ( F N ( r,z,ψ h ) /r ) ∂ψ is the derivative of F N ( r, z, ψ h ) /r with respect to its third argu-ment. Recall the definitions of the Gram matrix G in (3.5) and the test norm || · || V in (3.13), and we get G = (cid:18) G V G S (cid:19) , (3.17)where the block matrices are defined as( G V ) ij = (cid:16) e φ j , e φ i (cid:17) Ω h + (cid:16) (cid:101) ∇ · e φ j , (cid:101) ∇ · e φ i (cid:17) Ω h , ( G S ) ij = (cid:0) e τj , e τi (cid:1) Ω h + (cid:16) (cid:101) ∇ e τj , (cid:101) ∇ e τi (cid:17) Ω h . DAPTIVE DPG FOR GRAD-SHAFRANOV J T ( U ) G − (cid:104) B L U − B N ( U ) − F L (cid:105) = 0 . (3.18)The above system is the final formulation of the DPG scheme for the Grad-Shafranovequation. The system is nonlinear due to the source term and thus it requires specialcare to solve it efficiently. In Section 4 we will focus on exploring its nonlinear solvers.
4. Nonlinear solvers.
In this work, we have considered two types of nonlin-ear solvers, Newton’s method and Picard iteration equipped with Anderson acceler-ation, for the matrix-vector form (3.18). Some advantages and disadvantages arefound during the numerical experiments, which will be addressed in this section. Allthe nonlinear solvers discussed in this work are implemented through PETSc SNESsolvers [4].
The most popular method for nonlinear problems isNewton’s method, which is the first type of solvers we have investigated during thiswork. The challenge to implement a Newton’s method for the system (3.18) is toassemble its Jacobian matrix. It is easy to see that the Jacobian of the system (3.18)is given by J = J T ( U ) G − J ( U ) + H r . (4.1)Here r stands for the second component in the residual r ( u h ) and H is the Hessiandefined by − ∂D N ∂ Ψ where D N is the matrix associated with the nonlinear source term.To implement the Jacobian matrix J , one has to implement a Hessian-vector productassociated with H , which is not trivial especially in finite element packages such asMFEM. Instead, we avoid assembling the Hessian matrix by using the Jacobian-freeNewton-Krylov (JFNK) method that estimate J through a finite difference approxi-mation of J ≈ G ( U + (cid:15) V ) − G ( U ) (cid:15) , where the nonlinear function in this case is G ( U ) = J T ( U ) G − (cid:104) B L U − B N ( U ) − F L (cid:105) . Without a proper preconditioner, the JFNK solver is typically very inefficient.Through looking at its exact form of (4.1), it seems natural to provide a precondi-tioner of J T ( U ) G − J ( U ) for JFNK which only ignore one diagonal term related to H . A block Jacobi preconditioner is further provided to accelerate the inversion of J T ( U ) G − J ( U ) (the block Jacobi preconditioner is identical to the preconditionerfor J ( U ) G − B L that will be discussed in the next section and we therefore skip thediscussion here). A GMRES solver for the approximated J is further used in theouter iteration to guarantee the convergence of the full nonlinear solver.Such an approach works fine for problems that have very weak nonlinearity, whileit fails to converge for hard problems when the nonlinear source term is strong. Theignorance of the block matrix H r appears to be critical in those problems. Forinstance, for the nonlinear problems presented in Section 7, the JFNK solver with theabove preconditioning strategy fails, while the same solver works well for the linearproblems presented there. Therefore, we further explore another approach that isbased on Picard iteration, which is discussed in the next section.0 Z. PENG, Q. TANG, AND X.-Z. TANG
The second type of nonlinear solvers wehave investigated is the Picard iteration with the Anderson acceleration (Andersonmixing) method [2, 51]. We present the algorithm of Anderson acceleration that isused in this work
Algorithm 4.1
Anderson acceleration solving nonlinear fixed point problem U − A ( U ) = 0 . Given initial guess U . Set the initial residual R as R = U − A ( U ).Set U = (1 − λ ) U + λ A ( U ), R = U − A ( U ) and k = 1. while (cid:16) || r k || ≥ max( rtol || r || , atol ) and || x k − x k − || ≥ stol || x k − || (cid:17) do Set m k = min { k, m } .Find ( α k , . . . , α km k ) to be the solution of the constrained minimization problem:( α k , . . . , α km k ) = arg min (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) m k (cid:88) i =0 α ki R k − i (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) , such that m k (cid:88) i =0 α ki = 1 . Set U k +1 = (1 − λ k ) (cid:80) m k i =0 α ki U k − i + λ k (cid:80) m k i =0 α ki A ( U k − i ).Calculate the residual: R k +1 = U k +1 − A ( U k +1 ) and k := k + 1. end while Here λ k is determined by a cubic backtracking line search algorithm [19].Note that the system (3.18) is not in the form of a fixed point problem, and hencesome transformations are needed before applying the Anderson acceleration. Wehave investigated two types of rewritten systems. A straightforward rewritten formof (3.18) compatible with the Anderson acceleration leads to the following system U − (cid:110) U − J T ( U ) G − (cid:104) B L U − B N ( U ) − F L (cid:105)(cid:111) = 0 . However, we found that the nonlinear solver for this system converges poorly evenwith the Anderson acceleration. Instead, the following nonlinear fixed-point problemis used in our implementation U − (cid:0) J T ( U ) G − B L (cid:1) − (cid:104) J T ( U ) G − ( B N ( U ) − F L ) (cid:105) = 0 . (4.2)It is straightforward to see that (4.2) is equivalent with (3.18). Note that a similarstrategy was applied in [48, 49]. According to [5], this reformulation strategy can beviewed as a “preconditioning” procedure for the nonlinear solver.At each iteration, we need to invert the block matrix of J T ( U ) G − B L . Using(3.15), (3.16) and (3.17), with some work, we found that J T ( U ) G − B L = (cid:18) M r (cid:101) ∇ h T (cid:98) ψ div h − D N ( ψ h ) T (cid:99) q n (cid:19) T (cid:18) G V G S (cid:19) − (cid:18) M r (cid:101) ∇ h T (cid:98) ψ div h T (cid:99) q n (cid:19) = M Tr G − V M r + div Th G − S div h M Tr G − V (cid:101) ∇ h div Th G − S T (cid:99) q n M Tr G − V T (cid:98) ψ (cid:101) ∇ Th G − V M r − D N ( ψ h ) T G − S div h (cid:101) ∇ Th G − V (cid:101) ∇ h − D N ( ψ h ) T G − S T (cid:99) q n (cid:101) ∇ T G − V T (cid:98) ψ T T (cid:99) q n G − S div h T T (cid:99) q n G − S T (cid:99) q n T T (cid:98) ψ G − V M r T T (cid:98) ψ G − V (cid:101) ∇ h T T (cid:98) ψ G − V T (cid:98) ψ DAPTIVE DPG FOR GRAD-SHAFRANOV − D N ( ψ h ), and the flexible GMRES solver [46] is applied to invert this block matrix.Since the final form of J T ( U ) G − B L is rather complicated, it is necessary to providea few words on some details of matrix assembling when implementing the scheme.In our implementation, we assemble most of small blocks in the beginning of thesimulation except for two blocks associated with − D N ( ψ h ). During each iterationof the nonlinear solver, we only need to update two small block matrices resultingfrom the nonlinear source term. The inversion of the Gram matrices, G V and G S ,is implemented exactly by inverting the corresponding small matrix on each elementduring assembling. This is sufficient since our test space is discontinuous. The rest ofmatrix assembling uses standard approaches.The iterative linear solver to invert J T ( U ) G − B L is still inefficient without pre-conditioning. Since this inversion needs to be performed at each iteration of thenonlinear solver, it is required that such an iterative linear solver must be efficient.To improve the efficiency of this linear solver, we provide a block Jacobian precondi-tioner associated with J T ( U ) G − B L P = P P P P where the blocks are given by P := M Tr G − V M r + div Th G − S div h , P := (cid:101) ∇ Th G − V (cid:101) ∇ h ,P := T T (cid:99) q n G − S T (cid:99) q n , P := T T (cid:98) ψ G − V T (cid:98) ψ In the preconditioning stage, we inverted each block of P ii with algebraic multigridpreconditioners provided by HYPRE. We found it is sufficient to use the standardalgebraic multigrid method (AMG) [45] with one V-cycle to invert the blocks of P , P and P . The AMS method [31] with one V-cycle is applied to invert P . Notethat the block P is corresponding to (cid:98) q n,h , which is associated with the trace of a H (div) space. As an AMG method designed for H (div) and H (curl) spaces basedon the idea of auxiliary space method [25], the AMS method improves the efficiencyover the standard AMG method for those spaces. Hence, compared with the standardAMG method, we found that the AMS method is a better choice to invert the block P .
5. Adaptivity.
An adaptive mesh refinement (AMR) approach is developed inour implementation to improve the efficiency of our DPG solver. The basic idea ofAMR is to refine a mesh dynamically based on some estimate of local errors in thecurrent solution. For the steady state problem such as the case we consider in thiswork, AMR can improve the convergence of nonlinear solvers for some hard problemscompared with a uniformly refined mesh, since a significant amount of work is focusedon the region where errors are large in an adaptive mesh thanks to the error estimator.One advantage of the DPG method in the context of AMR is that it comes witha natural error estimator, while other numerical schemes typically rely on certainestimates, such as the Zienkiewicz-Zhu error estimator which estimates the errorsusing local numerical fluxes or error estimators through Richardson extrapolation. Forinstance, in [8] a natural posteriori error estimator based on the residual || B u − l || V (cid:48) is derived for the DPG method solving linear problems. In [7] the posteriori error2 Z. PENG, Q. TANG, AND X.-Z. TANG analysis is generalized to nonlinear problems and an AMR algorithm is presented forthe nonlinear problem.Here we first present the details of the error estimator for the DPG method andthen describe the details of an AMR algorithm based on it.
Recall (3.8), and the DPG method seeks u h such that u h = arg min w h ∈ U h ||B ω h − l || V (cid:48) . (5.1)A posteriori error estimate for DPG method solving linear problems is given in [8] as C ||B u h − l || V (cid:48) ≤ || u − u h || U ≤ C ||B u h − l || V (cid:48) + C osc( l ) , (5.2)where C , C and C are constants independent of the mesh. The term osc( l ) doesnot depend on the numerical solution u h , and it is a high order term with respectto the mesh size h , determined by the approximation properties of the finite elementspaces. With regularity assumptions, [7] generalizes the results to nonlinear problems.As shown in Section 3.1, ||B u h − l || V (cid:48) = sup v h ∈ V h , v h (cid:54) =0 | r ( u h , v h ) | || u h || V = r ( u h ) T G − r ( u h ) . (5.3)Using the Riesz representation, we can find ε u h ∈ V such that( ε u h , v ) V = b N ( u h , v ) − l ( v ) = r ( u h , v ) = ( B u h − l )( v ) , ∀ v ∈ V. Furthermore, || ε u h || V = ||B u h − l || V (cid:48) = (cid:113) r ( u h ) T G − r ( u h ). Then, we define the localerror estimator E K on each element K as E K := || ε u h || V ( K ) . (5.4)The corresponding total error estimator for the whole domain is defined as E total := (cid:115) (cid:88) K ∈ Ω h || ε u h || V ( K ) . (5.5) We only consider conforming h -adaptive refinementwith a fixed polynomial order in this work, although the non-conforming AMR isrecently developed in MFEM [10]. We use the following mesh-refinement strategyin the implementation. After each nonlinear solver for the DPG scheme, an updatedsolution u h for the Grad-Shafranov equation is obtained and its residual norm on eachelement can be estimated using (5.4). In our implementation, for a given element K ,it will be refined if all the following three conditions hold: (1) ε K > atol amr , (2) ε K > θ max max K (cid:48) ∈ Ω h ε K (cid:48) , (3) ε K > θ total E total / (cid:112) N mesh , where θ max ∈ [0 , θ total ∈ [0 ,
1) and atol amr ≥ K based on the criterion that its local error is lessthan a predetermined tolerance atol amr , and its estimated error is relatively largercompared to the errors in other elements. Finally, the stoppage criterion for the AMRiteration is when the total number of elements is larger than some upper bound, orno element satisfies the three conditions above. DAPTIVE DPG FOR GRAD-SHAFRANOV
6. Implementation.
The described DPG algorithm for the Grad-Shafranovequation is implemented under the framework of c++ library MFEM [1]. It would beappropriate to say a few words about our implementation under MFEM.The algorithm is implemented in parallel using a standard domain decompositionmethod. All the vectors and small block matrices use the parallel distributed datastructure provided by the package and its communication between sub-domains isbased on the message-passing interface (MPI). All the linear and nonlinear solversdescribed in this work are implemented through PETSc [4] and some of the smallmatrices are preconditioned with HYPRE [20] algebraic multigrid preconditioners toimprove efficiency. All the matrix assemblies are performed under MFEM to minimizethe communication cost between MFEM and PETSc.Our implementation of the DPG algorithm is general, supporting arbitrary or-der of accuracy and general meshes including triangular, quadrilateral and high-ordercurvilinear meshes, taking full advantage of the capabilities from the package. How-ever, our current implementation only supports conforming AMR, but the implemen-tation still offers refining and dynamic load-balancing, which will be performed afterthe mesh is updated in each AMR iteration. As a result, in Section 7, we mainly fo-cus on the numerical results on triangular meshes and its conforming adaptive meshrefinement to verify the scheme. Another reason to focus on triangular meshes is dueto the limitation of mesh generations. Currently, for all the problems with complexgeometry, we generate triangular meshes using the Gmsh software. To the best of ourknowledge, curvilinear meshes with complex geometry obtained from a mesh genera-tor are not fully supported by MFEM. Nevertheless, the focus of the current work ison the DPG scheme, AMR and efficient nonlinear solvers, while mesh generators arewell beyond the scope of this work. Therefore, it is sufficient to use triangular meshesgenerated by Gmsh for the numerical tests.
Reproducibility.
The implementation of the DPG scheme for the Grad-Shafranovequation, as well as the numerical examples presented in this work, is freely availableas an MFEM fork at github.com/ZhichaoPengMath/mfem .
7. Numerical tests.
We demonstrate the performance of our DPG schemethrough a set of numerical examples. Several examples are presented to study theaccuracy of the DPG scheme, including two linear examples and one nonlinear exam-ple. The linear examples are further used to compare the performance of the DPGscheme with other finite element methods. We then consider a nonlinear case in-volving slightly more complicated geometry. Finally, several nonlinear examples arepresented to study the performance of AMR.In practice, we use the trace of the Raviart-Thomas space [42, 37] for the spacecorresponding to the third component of U kh . For the discrete test space V k,sh , s = 2 ischosen throughout the numerical study. Unless otherwise noted, the results presentedhere use the Picard iteration with Anderson acceleration as the nonlinear solver. Two numerical examples are presented to demonstratethe accuracy of the DPG scheme. The first example considers the linear Solovevprofiles from [39, 9]. This test is further used to compare the accuracy for boththe solution ψ and its derivative (cid:101) ∇ ψ/r solved by our DPG method, conventionalcontinuous Galerkin (CG) method and the hybridized discontinuous Galerkin (HDG)method. The second example considers a manufactured solution with a nonlinearsource. We use it to demonstrate the accuracy of our scheme in nonlinear problems.4 Z. PENG, Q. TANG, AND X.-Z. TANG
The linear Grad-Shafranov equation of r (cid:101) ∇ · (cid:18) r (cid:101) ∇ ψ (cid:19) = r has an exact solution [39, 9] in the form of ψ ( r, z ) = r d + d r + d ( r − r z ) , where the parameters d , d and d determine the reasonable plasma cross section. G (1) of ITER G (1) of NSTX Fig. 1 . Computational domains and initial meshes of G (1) for two linear accuracy tests. Themeshes are generated through Gmsh by describing the computational boundaries. In [39] it has been described how to use the inverse aspect ratio ε , the elongation κ and the triangularity δ to determine the parameters d , d and d . We adopt thesame manufactured solutions in the accuracy test. The parameters are determinedby the following linear system ε ) (1 + ε ) − ε ) (1 − ε ) − δε ) (1 − δε ) − − δε ) κ ε d d d = − (1 + ε ) (1 − ε ) (1 − δε ) . Two tests are considered in this case, including the ITER-like [3] configuration of ε = 0 . κ = 1 . δ = 0 .
33 and NSTX-like configuration [47] of ε = 0 . k = 2, δ = 0 .
35. The computational domains and initial meshes are presented in Figure 1.The meshes are generated using Gmsh by describing the computational domains usingthe exact solutions. In particular, the computational boundaries are described byChebyshev nodes along the r direction and z coordinates are then determined bysolving ψ = 0. The numerical solutions are shown in Figure 2 and the presentedresults include ψ and its two derivatives of ψ r /r and ψ z /r . The results show a goodagreement with their exact solutions.This example is further used to study the accuracy of the DPG schemes as well asseveral commonly used finite element schemes. We presented the numerical solutionssolved by two other schemes, the HDG and CG methods. To approximate q = (cid:101) ∇ ψ/r , note that the ultraweak DPG and HDG methods directly solve q h in the weakformulations, while the CG method only solve ψ h and its q h is obtained by solving DAPTIVE DPG FOR GRAD-SHAFRANOV ITER ψ − . ITER ψ r /r − .
208 0 . ITER ψ z /r − .
120 0 . NSTX ψ − .
227 0
NSTX ψ r /r − .
542 0 . NSTX ψ z /r − .
269 0 . Fig. 2 . Numerical solutions and their derivatives for two linear accuracy tests in Section 7.1.1.The DPG scheme is used with the meshes of G (1) and quadratic polynomials. ( q h , φ ) = ( (cid:101) ∇ ψ h /r, φ ) for arbitrary φ in the discrete test space. It is well knownthat such an approach to compute q h will result in accuracy reduction by at leastone order. On the other hand, the results of the HDG method do not perform thepost-processing reconstruction, which typically could lift one order of accuracy for ψ . To demonstrate the results, as an example, quadratic polynomials are used forthe trial space of all three schemes. The uniform refinement is performed during theconvergence study and the initial meshes are given in Figure 1. L ∞ -errors of ψ and q are presented for different schemes in Table 1. Here the L ∞ -error is defined as L ∞ error ( u h ) = max K ∈ Ω h || u h ( x ) − u exact ( x ) || ∞ ,K , where ||·|| ∞ ,K is the standard L ∞ norm on K . The L ∞ -error for the vector q is takenas the maximal error among all its components. All three schemes achieve a third-order accuracy for ψ as expected. The accuracy for q of the DPG and HDG methodsare third order, while the CG method only has second order accuracy. Compared withthe HDG method, though with the same accuracy order, the error of the DPG methodis smaller for both ψ and q . We note that the HDG method is computationally moreefficient than the DPG scheme for this case, as it results into a linear system of smallersize through a Schur complement. The Schur complement for the matrix of B T G − B in the linear DPG scheme is not as obvious as the HDG scheme. It becomes evenmore challenging for the nonlinear problems we will consider later. Therefore, we donot investigate this direction in this work. It is also important to mention that thepost-processing reconstruction for the HDG scheme can not improve the accuracy of q = (cid:101) ∇ ψ/r , though it lifts the accuracy order of ψ .6 Z. PENG, Q. TANG, AND X.-Z. TANG
The convergence results indicate the DPG scheme produce more accurate derivatesthan the other two schemes.
Table 1 L ∞ -errors and orders of ψ and q = − (cid:101) ∇ ψ/r for the linear problems in Section 7.1.1. TheDPG, HDG and CG schemes are used with P polynomials. It demonstrates the accuracy of theDPG scheme and also shows the performance of the different schemes when computing its derivative q . L ∞ -errors and orders of ψ Test Grid DPG HDG CG ψ order ψ order ψ order G (1) G (2) G (4) G (8) G (1) G (2) G (4) G (8) L ∞ -errors and orders of q Test Grid DPG HDG CG q order q order q order G (1) G (2) G (4) G (8) G (1) G (2) G (4) G (8) Next consider a manufactured solution given in [48].The solution satisfies ψ ( r, z ) = sin ( k r ( r + r )) cos( k z z ) , with the source term given by F ( r, z, ψ ) =( k r + k z ) ψ + k r r cos( k r ( r + r )) cos( k z z ) + r (cid:104) sin ( k r ( r + r )) cos ( k z z ) − ψ + exp (cid:0) − sin (cid:0) k r ( r + r ) (cid:1) cos( k z z ) (cid:1) − exp( − ψ ) (cid:105) , where the coefficients are k r = 1 . π , k z = 1 .
15 and r = − .
5. We consider the sameITER geometry as the first linear test in the previous section. The computationaldomain and the initial mesh are identical to the ITER case given in Figure 1. For thiscase, a non-homogenous Dirichlet boundary condition given by the exact solution isused throughout the computations. The solutions are presented in Figure 3, in whichquadratic polynomials are used for the trial space.
DAPTIVE DPG FOR GRAD-SHAFRANOV U kh with k = 1 , , L ∞ numerical errors and the corresponding orders are presentedin Table 2. For the discrete trial space U kh with k = 1 ,
2, the optimal convergenceof ( k + 1)-th order is observed for both the solutions and their derivatives. For U h ,though order reduction of ψ occurs, optimal convergence for its derivatives is stillobserved. ψ .
234 1 ψ r /r − .
23 4 . ψ z /r − .
646 0 . Fig. 3 . Numerical solutions and their derivatives for the nonlinear accuracy problem in Sec-tion 7.1.2. The DPG scheme is used with the meshes of G (1) and quadratic polynomials. Table 2 L ∞ errors and orders of ψ and q = − (cid:101) ∇ ψ/r for the nonlinear problem in Section 7.1.2. P , P and P polynomial spaces are used for the DPG scheme. L ∞ -errors and orders of ψ Grid P P P G (1) G (2) G (4) G (8) L ∞ -errors and orders of q Grid P P P G (1) G (2) G (4) G (8) Next we consider a D-shaped geometry taken from [48].The boundary of computational domain ω is determined by r ( s ) = 1 + ε cos( s + arcsin( δ sin( s )) , z ( s ) = εκ sin( s ) , Z. PENG, Q. TANG, AND X.-Z. TANG where s ∈ [0 , π ], ε = 0 . δ = 0 .
33 and κ = 1 .
7. The source term is F ( r, z, ψ ) = r (cid:104) −
12 (1 − ψ ) (cid:105) . (7.1)A homogenous Dirichlet boundary condition is used in the simulation. Numerical re-sults with trial space U h are shown in Figure 4. This case demonstrates the capabilityto handle different geometry and nonlinear source in our DPG implementation. ψ . . ψ r /r − . . ψ z /r − . . Fig. 4 . Numerical solutions for the case of the D-shaped geometry in Section 7.2.
In the final part of the numerical example, we focus on theperformance of the DPG scheme when the AMR strategy is applied. Two tests arepresented with one involving a rectangular geometry and one involving D-shapedgeometry. We set θ max = 0 .
025 and θ global = 0 . amr and thediscrete trial space U kh will be specified in each test. Initial mesh Mesh after 4 adaptive refinement Errors vs element number
Fig. 5 . Left: initial mesh for the test in Section 7.3.1. Middle: mesh after four adaptiverefinement iterations. Right: errors verse total number of elements for both uniform and adaptiverefinements.
First consider a rectangle computational do-main of Ω = [0 . , . × [ − . , .
75] with a homogenous boundary condition of ψ = 0 .
25. The source term is nonlinear given by F ( r, ψ ) = 2 r ψ (cid:104) c (1 − exp( − ψ /σ ) + 1 σ ( c + c ψ ) exp( − ψ /σ ) (cid:105) , DAPTIVE DPG FOR GRAD-SHAFRANOV σ = 0 . c = 0 . c = 0 .
2. The reason for first considering therectangular geometry is to eliminate the impact of the boundary condition and theextra complexity associated with a complex geometry when performing AMR. Thefocus of this case is therefore on the improvement of numerical errors of AMR over auniform mesh refinement.For this test, we set atol amr =1e-8 and use quadratic polynomial space U h asthe discrete trial space. The initial mesh is presented in Figure 5. Note that theinitial mesh is chosen carefully so that it retains a symmetry along y = 0, which isconsistent with the exact solution. In our solver, we perform one refinement checkafter each nonlinear solve, and the updated adaptive mesh will be used again in thenext nonlinear solve. An iteration procedure (which is referred to as AMR iterationsin this work) is performed until a targeting error goal is achieved or the total numberof elements is larger than an upper bound.The mesh after four AMR iterations is presented in Figure 5. It is accompaniedby the numerical solutions on the initial and adaptive meshes presented in Figure 6.The color in the adaptive mesh of Figure 5 indicates different refinement levels ofAMR. We first find that the adaptive mesh remains symmetric throughout the AMRiterations. Note that q = (cid:101) ∇ ψ/r has sharp features near the top-right corner and thebottom-right corner, which are not well-resolved on the initial mesh. After four AMRiterations, these features are better captured and most refinements are performed nearthem. Overall, the computational efforts are focused on the region where interestingphysics happens during the AMR iteration. All those facts indicate our AMR strategyis effective and efficient.As a further verification of the AMR strategy, we also compare the numericalerrors under uniform and adaptive refinements. Here the numerical error is measuredby the value of our error estimator defined in (5.5). The convergence histories forboth AMR and uniform refinement are also presented in Figure 5. The numericalerrors are presented under increasing total number of elements. Note that the errorson AMR are much smaller than errors on uniform meshes when the total number ofelements (degree of freedom) are the same. Another important observations is thatas the meshes are refined, the numerical errors on AMR decay quicker than the errorson uniform meshes. Therefore, the convergence histories quantitatively confirm theefficiency of the AMR approach in this work. To further confirm the effectiveness of AMR onproblems involving complex geometry, we consider the same nonlinear test in D-shaped geometry as described in Section 7.2. For this test, we set atol amr =1e-6 anduse linear polynomial space U h as the discrete trial space for the DPG scheme.The initial mesh and the convergence history of AMR and uniform refinement arepresented in Figure 7. We observe that, compared with the uniform refinement, thenumerical error on AMR is much smaller and decays faster. This example implies theefficiency and effectiveness of our AMR strategy on problems with complex geometry.
8. Conclusions.
This work focuses on designing and developing an arbitrary-order adaptive DPG method for the nonlinear Grad-Shafranov equation. The mainfocus of the work is to investigate advantages of the DPG scheme when solving thederivatives as the magnetic field in the Grad-Shafranov equation. The ultraweakformulation of the DPG scheme is developed based on a minimal residual method andthe efficient nonlinear solvers are proposed and studied. The algorithm is augmentedwith an AMR strategy to improve the efficiency of the algorithm. The proposedalgorithm is implemented in parallel under the framework of MFEM.0
Z. PENG, Q. TANG, AND X.-Z. TANG
Initial mesh ψ .
249 0 . Initial mesh ψ r /r − . . Initial mesh ψ z /r − . . After 4 refinement ψ .
251 0 . After 4 refinement ψ r /r − . . After 4 refinement ψ z /r − . . Fig. 6 . Numerical solutions and their derivatives for the rectangular geometry case in Section7.3.1. The first row shows the results on the initial mesh. The second row shows the results afterfour AMR iterations.
Initial mesh Errors vs element number
Fig. 7 . Left: initial mesh for the test in Section 7.3.2. Right: errors verse total number ofelements for both uniform and adaptive refinements.
A series of numerical results is presented to verify the accuracy and efficiency ofthe algorithm. In particular, the algorithm is found to produce more accurate numer-ical derivatives than some of the commonly used finite element schemes such as CGand HDG schemes. It is also demonstrated that the scheme produces optimal con-vergence for both the solution and derivatives for k -th order polynomial spaces. Thenonlinear solver based on Anderson iteration is found to be efficient to solve problemsinvolving strong nonlinearity. Finally, numerical examples also quantitatively demon-strate the improvement of efficiency through the AMR strategy compared with the DAPTIVE DPG FOR GRAD-SHAFRANOV
REFERENCES[1]
MFEM: Modular finite element methods library . .[2]
D. G. Anderson , Iterative procedures for nonlinear integral equations , Journal of the ACM(JACM), 12 (1965), pp. 547–560.[3]
R. Aymar, P. Barabaschi, and Y. Shimomura , The ITER design , Plasma physics and con-trolled fusion, 44 (2002), p. 519.[4]
S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin,A. Dener, V. Eijkhout, W. Gropp, et al. , PETSc users manual , (2019).[5]
P. R. Brune, M. G. Knepley, B. F. Smith, and X. Tu , Composing scalable nonlinear alge-braic solvers , SIAM Review, 57 (2015), pp. 535–565.[6]
T. Bui-Thanh and O. Ghattas , A PDE-constrained optimization approach to the discontin-uous Petrov–Galerkin method with a trust region inexact Newton-CG solver , ComputerMethods in Applied Mechanics and Engineering, 278 (2014), pp. 20–40.[7]
C. Carstensen, P. Bringmann, F. Hellwig, and P. Wriggers , Nonlinear discontinuousPetrov–Galerkin methods , Numerische Mathematik, 139 (2018), pp. 529–561.[8]
C. Carstensen, L. Demkowicz, and J. Gopalakrishnan , A posteriori error control for DPGmethods , SIAM Journal on Numerical Analysis, 52 (2014), pp. 1335–1353.[9]
A. J. Cerfon and J. P. Freidberg , One size fits all analytic solutions to the Grad–Shafranovequation , Physics of Plasmas, 17 (2010), p. 032502.[10]
J. ˇCerven`y, V. Dobrev, and T. Kolev , Non-conforming mesh refinement for high-orderfinite elements , arXiv preprint arXiv:1905.04033, (2019).[11]
J. Chan, L. Demkowicz, and R. Moser , A DPG method for steady viscous compressible flow ,Computers & Fluids, 98 (2014), pp. 69–90.[12]
J. Chan, L. Demkowicz, R. Moser, and N. Roberts , A class of discontinuous petrov–galerkinmethods. part v: Solution of 1d burgers and navier–stokes equations , ICES Report, 29(2010).[13]
L. Demkowicz and J. Gopalakrishnan , A class of discontinuous petrov–galerkin methods.part i: The transport equation , Computer Methods in Applied Mechanics and Engineering,199 (2010), pp. 1558–1572.[14]
L. Demkowicz and J. Gopalakrishnan , A class of discontinuous Petrov–Galerkin methods.II. Optimal test functions , Numerical Methods for Partial Differential Equations, 27 (2011),pp. 70–105.[15]
L. Demkowicz and J. Gopalakrishnan , Analysis of the dpg method for the poisson equation ,SIAM Journal on Numerical Analysis, 49 (2011), pp. 1788–1809.[16]
L. Demkowicz and J. Gopalakrishnan , Discontinuous Petrov–Galerkin (DPG) Method , En-cyclopedia of Computational Mechanics Second Edition, (2017), pp. 1–15.[17]
L. Demkowicz, J. Gopalakrishnan, and A. H. Niemi , A class of discontinuous Petrov–Galerkin methods. Part III: Adaptivity , Applied numerical mathematics, 62 (2012),pp. 396–427.[18]
L. F. Demkowicz and J. Gopalakrishnan , An overview of the discontinuous Petrov Galerkinmethod , in Recent developments in discontinuous Galerkin finite element methods for par-tial differential equations, Springer, 2014, pp. 149–180.[19]
J. E. Dennis Jr and R. B. Schnabel , Numerical methods for unconstrained optimization andnonlinear equations , vol. 16, Siam, 1996.[20]
R. D. Falgout and U. M. Yang , hypre: A library of high performance preconditioners , inInternational Conference on Computational Science, Springer, 2002, pp. 632–641.[21] J. Gopalakrishnan and W. Qiu , An analysis of the practical DPG method , Mathematics ofComputation, 83 (2014), pp. 537–552.[22]
H. Grad and H. Rubin , Hydromagnetic equilibria and force-free fields , Journal of NuclearEnergy (1954), 7 (1958), pp. 284–285.[23]
Z. Guo, X.-Z. Tang, and C. J. McDevitt , Models of primary runaway electron distributionin the runaway vortex regime , Physics of Plasmas, 24 (2017), p. 112508.[24]
H. Heumann and F. Rapetti , A finite element method with overlapping meshes for free-boundary axisymmetric plasma equilibria in realistic geometries , Journal of ComputationalPhysics, 334 (2017), pp. 522–540. Z. PENG, Q. TANG, AND X.-Z. TANG[25]
R. Hiptmair and J. Xu , Nodal auxiliary space preconditioning in H (curl) and H (div) spaces ,SIAM Journal on Numerical Analysis, 45 (2007), pp. 2483–2509.[26]
E. Howell and C. R. Sovinec , Solving the grad–shafranov equation with spectral elements ,Computer Physics Communications, 185 (2014), pp. 1415–1421.[27]
G. Huysmans and O. Czarny , MHD stability in X-point geometry: simulation of ELMs ,Nuclear fusion, 47 (2007), p. 659.[28]
G. Huysmans, J. Goedbloed, W. Kerner, et al. , Isoparametric bicubic hermite elementsfor solution of the grad-shafranov equation , International Journal of Modern Physics C, 2(1991), pp. 371–376.[29]
S. Jardin , Computational methods in plasma physics , CRC Press, 2010.[30]
B. Keith, S. Petrides, F. Fuentes, and L. Demkowicz , Discrete least-squares finite elementmethods , Computer Methods in Applied Mechanics and Engineering, 327 (2017), pp. 226–255.[31]
T. V. Kolev and P. S. Vassilevski , Parallel auxiliary space AMG for H (curl) problems ,Journal of Computational Mathematics, (2009), pp. 604–623.[32]
J. Lee and A. Cerfon , Ecom: A fast and accurate solver for toroidal axisymmetric mhdequilibria , Computer Physics Communications, 190 (2015), pp. 72–88.[33]
H. L¨utjens, A. Bondeson, and O. Sauter , The chease code for toroidal mhd equilibria ,Computer physics communications, 97 (1996), pp. 219–260.[34]
D. Moro, N. Nguyen, and J. Peraire , A hybridized discontinuous Petrov–Galerkin schemefor scalar conservation laws , International Journal for Numerical Methods in Engineering,91 (2012), pp. 950–970.[35]
S. Nagaraj, J. Grosek, S. Petrides, L. F. Demkowicz, and J. Mora , A 3D DPG Maxwellapproach to nonlinear Raman gain in fiber laser amplifiers , Journal of ComputationalPhysics: X, 2 (2019), p. 100002.[36]
S. Nagaraj, S. Petrides, and L. F. Demkowicz , Construction of DPG Fortin operators forsecond order problems , Computers & Mathematics with Applications, 74 (2017), pp. 1964–1980.[37]
J.-C. N´ed´elec , Mixed finite elements in R , Numerische Mathematik, 35 (1980), pp. 315–341.[38] A. Palha, B. Koren, and F. Felici , A mimetic spectral element solver for the grad–shafranovequation , Journal of Computational Physics, 316 (2016), pp. 63–93.[39]
A. Pataki, A. J. Cerfon, J. P. Freidberg, L. Greengard, and M. ONeil , A fast, high-ordersolver for the Grad–Shafranov equation , Journal of Computational Physics, 243 (2013),pp. 28–45.[40]
S. Petrides , Adaptive multilevel solvers for the discontinuous Petrov–Galerkin method withan emphasis on high-frequency wave propagation problems , PhD thesis, 2019.[41]
S. Petrides and L. F. Demkowicz , An adaptive DPG method for high frequency time-harmonic wave propagation problems , Computers & Mathematics with Applications, 74(2017), pp. 1999–2017.[42]
P.-A. Raviart and J.-M. Thomas , A mixed finite element method for 2nd order elliptic prob-lems , in Mathematical aspects of finite element methods, Springer, 1977, pp. 292–315.[43]
N. V. Roberts, T. Bui-Thanh, and L. Demkowicz , The DPG method for the Stokes problem ,Computers & Mathematics with Applications, 67 (2014), pp. 966–995.[44]
N. V. Roberts, L. Demkowicz, and R. Moser , A discontinuous Petrov–Galerkin methodologyfor adaptive solutions to the incompressible Navier–Stokes equations , Journal of Compu-tational Physics, 301 (2015), pp. 456–483.[45]
J. W. Ruge and K. St¨uben , Algebraic multigrid , in Multigrid methods, SIAM, 1987, pp. 73–130.[46]
Y. Saad , A flexible inner-outer preconditioned GMRES algorithm , SIAM Journal on ScientificComputing, 14 (1993), pp. 461–469.[47]
S. Sabbagh, S. Kaye, J. Menard, F. Paoletti, M. Bell, R. Bell, J. Bialek, M. Bitter,E. Fredrickson, D. Gates, et al. , Equilibrium properties of spherical torus plasmas inNSTX , Nuclear Fusion, 41 (2001), p. 1601.[48]
T. S´anchez-Vizuet and M. E. Solano , A Hybridizable Discontinuous Galerkin solver for theGrad–Shafranov equation , Computer Physics Communications, 235 (2019), pp. 120–132.[49]
T. S´anchez-Vizuet, M. E. Solano, and A. J. Cerfon , Adaptive Hybridizable Discontinu-ous Galerkin discretization of the Grad-Shafranov equation by extension from polygonalsubdomains , arXiv preprint arXiv:1903.01724, (2019).[50]
V. Shafranov , On magnetohydrodynamical equilibrium configurations , Soviet Physics JETP,6 (1958), p. 1013.[51]