An Anisotropic p-Adaptation Multigrid Scheme for Discontinuous Galerkin Methods

In this work, we present an anisotropic p-adaptation multigrid algorithm for discontinuous Galerkin methods for steady-state problems that uses a p-multigrid scheme both as a solver and as an anisotropic error estimator. To achieve this, we develop a new anisotropic truncation error estimator based on the tau-estimation method, that can be evaluated inside the multigrid cycle with a negligible extra cost. The new error estimator is cheaper to evaluate and more accurate than previous versions of the tau-estimation procedure. In our technique, a non-converged solution in a reference mesh is used to estimate the truncation error with the multigrid scheme for different combinations of polynomial orders in different directions inside every element, and the mesh is adapted accordingly to target a desired truncation error threshold. The accuracy and computational cost of the proposed p-anisotropic adaptation algorithm is tested for the steady viscous flow past a NACA0012 airfoil. A speed-up of 16 can be achieved in the proposed numerical example compared with the uniformly refined simulation without multigrid.

Multigrid methods speed up the iterative solution of large systems of equations using coarse-grid representations (lower levels). Iterative methods (known as smoothers in the multigrid community) are good at eliminating the high frequencies of the error fast; therefore, when applied to coarse-grid representations, they also reduce the low frequencies of the error. They have been broadly used in the high-order community in recent years in the form of p-multigrid [8,5] (where levels are constructed using different polynomial orders) and hp-multigrid [21,14] (where both the order and size of the elements are changed). Two types of multigrid methods can be found in the literature: linear and nonlinear multigrid. In our work, we make use of the nonlinear multigrid scheme, also known as the Full Approximation Scheme (FAS), since it enables the estimation of the truncation error of coarse representations, as will be shown. The smoother can be either a time-marching scheme (implicit or explicit), or an iterative method applied to the linearized problem.
Because of the allowed discontinuities on element interfaces, DG methods are capable of handling non-conforming meshes with hanging nodes and/or different polynomial orders efficiently [15,13,7]. It is possible to take advantage of this feature to accelerate the computations through local adaptation strategies. Local adaptation can be performed by subdividing or merging elements (h-adaptation) or by enriching or reducing the polynomial order in certain elements (p-adaptation). The main idea behind these methodologies is to reduce the number of degrees of freedom (NDOF) while maintaining a high accuracy, which translates into shorter computational times and reduced storage requirements. Furthermore, since several 2D and 3D implementations of the DG methods use tensor-product basis functions, it is possible to adapt the polynomial order in each coordinate direction independently. In order to identify the localized regions that need increased or decreased accuracy, an error estimator is commonly used.
There are several approaches to estimate the error and drive an adaptation method. In this work, we focus on truncation error estimates since it has been shown that a reduction of the truncation error controls the numerical accuracy of all functionals [10], hence reducing the truncation error necessarily leads to a more accurate lift and drag. The τ-estimation method [2] is a way to estimate the truncation error locally that has been used to drive mesh adaptation strategies in low-order [9,20] and high-order methods [17, 10,18]. The adaptation strategy consists in converging a high order representation (reference mesh) to a specified global residual and then performing a single error estimation followed by a corresponding mesh adaptation process. Rueda-Ramírez et al.
[19] developed a new method for estimating the truncation error of anisotropic representations that is cheaper to evaluate than previous implementations, and showed that it produces very accurate extrapolations of the truncation error, which enables the use of coarser reference meshes.
In this work, we employ the anisotropic truncation error estimator developed in [19] and the anisotropic p-adaptation method detailed in [18] to accelerate the computation of the compressible steady viscous flow past a NACA0012 at angle of attack 5 o , Re ∞ = 200 based on the airfoil chord, and M ∞ = 0.2. This particular settings correspond to a steady laminar flow, but the proposed method can be directly used with any steady solution (e.g. RANS). The paper is organized as follows: In section 2, we briefly describe the methods used in this paper. In section 3, we compare the performance of the proposed methods with traditional strategies for solving the flow past a NACA0012 and show the speed-up advantages for different accuracies. Finally, the conclusions are summarized in section 4.

DG Method
We consider the approximation of systems of conservation laws, where q is the vector of conserved variables, F is the flux dyadic tensor, and s is a source term. The domain Ω is partitioned in a mesh T = {e} consisting of K non-overlapping elements Ω e . Multiplying equation (1) by a test function v and integrating by parts over each subdomain Ω e yields the weak formulation: Let q, s, F and v be approximated by piece-wise polynomial functions defined in the space of L 2 functions: is the space of polynomials of degree at most N. The functions in V N can be represented in each element as a linear combination of basis functions where φ N i are usually tensor product expansions. After some manipulations, the discontinuous Galerkin finite element discretization system is obtained: where [M] is the mass matrix and F is a nonlinear operator, which are the assembled global versions of the element-wise mass matrices and nonlinear operators: where F F F e i is the i th position of the vector F F F e , which contains the value of F e for all the degrees of freedom of element e. In the rest of this paper, bold uppercase Roman letters and bold Greek letters are used to note vectors spanning several degrees of freedom, unless specified.
The numerical flux function F * allows to uniquely define the flux at the element interfaces and to weakly prescribe the boundary data as a function of the conserved variable on both sides of the boundary/interface and the normal vector. In the present work, we use Roe [16] as the advective Riemann solver and Bassi-Rebay 1 [4] as the diffusive Riemann solver.

Full Approximation Scheme p-Multigrid
The Full Approximation Scheme (FAS) is a nonlinear version of the multigrid method that is specially suited to solve systems of nonlinear equations [2]. Departing from equation (3) and defining the operator After β 1 sweeps of a smoother, a non-converged solutionQ P is obtained that has an associated discretization error P = Q P −Q P . The FAS multigrid procedure consists in obtaining an approximation to the discretization error in a coarse grid of order N and projecting it to the original problem of order P: where I P N is an L 2 projection operator N → P and Q N is the solution to the coarse-grid problem: where the source term is defined as In practice, several p-multigrid levels are used in V-or W-cycles. The smoothing steps that are performed when coarsening are called pre-smoothing sweeps, and the ones performed when refining back are called post-smoothing sweeps. Furthermore, Q N is not obtained exactly in the coarse grids, but approximated using an iterative methodQ N → Q N . In this work, we use a third order low-storage Runge-Kutta (RK3) as the smoother and V-cycles.

τ-Based p-adaptation
In this section we show how to drive an anisotropic p-adaptation procedure using the truncation error, which is estimated in the multigrid procedure.

The Anisotropic τ-Estimation Method
The non-isolated truncation error of a discretization of order N is defined as where q is the exact solution to the problem, I N is a discretizing operator, R is the continuous partial differentiation operator, and R N is the discrete partial differentiation operator. From equations (1) and (3): where I I I N is an operator that samples the exact solution on the points that correspond to the degrees of freedom of a representation of order N, and therefore equation (12) corresponds to the sampled values of R N (I N q).
Note that in steady cases, R(q) = 0 holds. Since the exact solution q is usually not at hand, we utilize the quasi a-piori τ-estimation method, which approximates the exact solution with the non-converged solution on a high-order grid q ≈q P , where N < P. Therefore, the steady non-isolated truncation error estimation yields On the left side of the arrow is the estimation of the truncation error that lives in the space V N , and on the right side is the sampled form of the truncation error estimation on the points that correspond to the degrees of freedom. In a DG representation, one can also define the isolated truncation errorτ aŝ whereF is the assembled version of the isolated nonlinear operator, defined elementwise as Note that equation (15) is (5) without substituting F by the numerical flux F * . This change eliminates the influence of the neighboring elements and boundaries on the truncation error of each element. We drop the hat notation in the next statements since they are valid for both the isolated and non-isolated truncation error.
The τ-estimation method can also be used with anisotropic representations, i.e.
where N i and P i are the polynomial orders in the direction i of the analyzed representation and the high-order reference solution, respectively, where N i < P i . Addi- Fig. 1 (a) Truncation error map for a specific element that shows log τ N1 N2 5,5 ∞ as a function of N 1 and N 2 (the black box shows the limit between the estimated and extrapolated maps). (b) Map of degrees of freedom (the black boxes show the polynomial orders that achieve τ N1 N2 tionally, Rueda-Ramírez et al. [19] showed that the truncation error of an anisotropic representation can be estimated using directional components: where the directional components in discrete form are therefore, and that these directional components decrease exponentially with the polynomial order in smooth solutions. Consequently, it is possible to use a semi-converged solutionq P 1 P 2 to estimate τ N 1 N 2 (N i < P i ) and then extrapolate the directional components τ i to obtain the values of τ N 1 N 2 for N i > P i . Fig 1(a) shows a graphical representation of the truncation error τ N 1 N 2 as estimated with a semi-converged solution of order P 1 = P 2 = 5.

The p-Adaptation Multigrid Scheme
It has been shown that the use of FAS p-multigrid methods speeds up the computation of steady-state and unsteady solutions of the compressible Navier-Stokes equations [5,8]. In addition, Rueda-Ramírez et al. [18] showed that the truncation error of an anisotropic representation can be inexpensively obtained inside an anisotropic p-multigrid cycle that performs the coarsening in one coordinate direction at a time.
In fact, the second term of equation (18) is naturally computed in an anisotropic multigrid for obtaining the coarse-grid source term (equation (9)). Therefore, we propose a p-adaptation multigrid scheme that makes use of the multigrid as a solver, but also as an error estimator. Every time the error is estimated, an anisotropic p-multigrid strategy is used to generate a truncation error map for each element, like the one in Fig 1(a). Afterwards, the polynomial orders in the different coordinate directions are selected for each element, such that a truncation error threshold τ max is achieved with the minimum NDOF possible, as illustrated in Fig. 1(b). In the simulations shown in this paper, the reference representation,q P , is converged to a residual τ max /10 before the p-adaptation stage, so that the truncation error is accurately estimated down to τ max , as was shown necessary by Kompenhans et al. [10].

Flow Past a NACA0012 Airfoil
In this section, we compare the performance of the proposed p-adaptation multigrid scheme with a uniformly adapted p-multigrid method (without local p-adaptation) and a uniformly adapted RK3 method when solving the steady viscous flow past a NACA0012 airfoil at angle of attack 5 o , Re ∞ = 200 (L ∞ = L chord ) and M ∞ = 0.2. This particular settings correspond to a steady laminar flow, but the proposed method can be directly used with any steady solution (e.g. RANS). An unstructured mesh of 2011 quadrilateral elements is employed (Fig. 2).
In the cases where multigrid is employed, the RK3 scheme is used as the iterative method (smoother), so that additional speed-ups are only due to the methods exposed in section 2. As in [18], a residual-based smoothing strategy is performed. The minimum number of smoothing sweeps is β = 200 for the coarsest multigrid level (N = 1) and β = 50 for any other level. After every β pre-smoothing sweeps, the residual in the next (coarser) representation is checked. If R N ∞ < 1.2 R N −1 ∞ , the pre-smoothing is stopped; otherwise, β additional sweeps are performed. Similarly, the norm of the residual after the post-smoothing is forced to be at least as low as it was after the pre-smoothing, R N post ∞ ≤ R N pre ∞ . If that condition is not fulfilled, additional β sweeps are taken until it is.
The isolated truncation error estimate is used to drive the p-adaptation method since it has been shown to provide better results than the non-isolated one [17,18,19]. The conservative form (equation (1)) of the compressible Navier-Stokes equations is discretized using the Discontinuous Galerkin Spectral Element Method (DGSEM) [1,11], which is a nodal (collocation) version of a DG method that uses Gauss points as the solution nodes and quadrature points, obtaining diagonal mass matrices. However, the methods that are exposed here can be applied to any DG scheme with tensor-product basis functions.
In [18] it was explained that, when using the DGSEM in general 3D curved meshes and p-nonconforming representations, the order of the mapping must be at most M ≤ N/2 for the numerical representation to be free-stream preserving. For this reason, the use of a conforming algorithm was proposed, which forces the polynomial orders to be conforming in the first layer of elements on a curved boundary. The use of a conforming algorithm is necessary to retain the well-known M ≤ N condition of the DGSEM [12]. In this work, we use the conforming algorithm on the airfoil surface since it showed to produce better results, although its use is not imperative as the considered test case is 2D.
For the uniformly adapted cases, the polynomial order is varied between N = 2 and N = 7. For the cases with local p-adaptation, a single-stage anisotropic p-adaptation procedure is performed, and the minimum polynomial order after adaptation is set to  N min = 1, whereas the maximum polynomial order after adaptation is set to N max = 7. The relative drag and lift errors of the adapted meshes are assessed by comparing with a reference solution of order N = 8: Fig. 3 shows a comparison between the errors obtained using theτ-based adaptation procedure and the ones using uniform p-refinement. As can be observed, the number of degrees of freedom is substantially reduced for the same accuracy when using theτ-based p-adaptation. This reduction translates into a reduction of the CPU-times. It is interesting to point out that, as the isolated truncation error thresholdτ max is decreased, the polynomial orders of the mesh tend to the maximum specified polynomial order, N max = 7. Consequently, the lift and drag coefficients also tends to C N =7 l . Using Fig. 3, it is possible to compute a speed-up for different levels of accuracy. Table 1 summarizes the speed-up calculations for the maximum level of accuracy that was achieved for the drag and lift coefficients.  Fig. 4 shows the distribution of polynomial orders after the single-stage adaptation procedure for a threshold of τ max = 5 × 10 −4 , which has related errors of e N =8 drag = 4.10 × 10 −5 and e N =8 lift = 7.31 × 10 −5 . As can be observed, the elements that are  enriched are mainly the ones on the boundary layer (specially leading and trailing edge), and the zones of the wake where the element size changes significantly.

Conclusions
In this work, we have applied recently developed error estimators and anisotropic p-adaptation methods in conjunction with multigrid solving strategies for solving the compressible Navier-Stokes equations. In particular, we have shown that the coupling of anisotropic truncation error-based p-adaptation methods with p-multigrid schemes can speed up the computation of steady-state solutions of PDEs. The achieved speedup depends on the desired accuracy, being this method optimal when high accuracy is required (low errors). In particular, a speed-up of 16.13 was achieved for the computation of the steady compressible viscous flow past a NACA0012 airfoil at angle of attack 5 • with respect to the uniformly adapted representation without multigrid.