Continuous adjoint-based error estimation and its application to adaptive discontinuous Galerkin method

An adaptive mesh refinement algorithm based on a continuous adjoint approach is developed. Both the primal equation and the adjoint equation are approximated with the discontinuous Galerkin (DG) method. The proposed adaptive algorithm is used in compressible Euler equations. Numerical tests are made to show the superiority of the proposed adaptive algorithm.


Introduction
In lots of applications, e.g., engineering analysis and design, engineers are interested in not only the whole solution but also its certain values, i.e., output functionals. For example, in computational aerodynamics, lift and drag are usually more important than the whole solution in the entire space [1] . In such computations, engineers are interested in the accuracy of these output functionals. Considering the finite computational resources, the adaptive approach is an efficient way to improve the calculation precision of these functionals. The adaptive approach usually uses local error indicators to guide the mesh refinement. The indicators derived based on the solution features such as the large gradient in flow variables and flow discontinuities have been widely used in adaptive algorithms. Although the feature-based adaptation has had some successes, it is liable to overrefine the regions that have little effect on the output functionals [2][3] .
In order to derive reliable and efficient adaptive algorithms to produce the desired increase in the accuracy of the output functionals, adjoint-based adaptive methods have been widely developed and studied [4][5][6][7][8][9][10][11][12][13] . Becker and Rannacher [4][5] proposed the dual weighted residual method based on the property of Galerkin orthogonality of the finite element method. Venditti [6] and Venditti and Darmofal [7] derived the adaptive procedure by the local error indicators constructed from both the primal solution and the adjoint solution in the standard finite volume discretization. Hartmann and Honston [8][9] and Hartmann [10] made extensive studies on the adaptive discontinuous Galerkin (DG) method based on an adjoint approach for aerodynamic flows. Wang and Mavriplis [11] developed the h-p adaptation approach for compressible Euler equations. Fidkowski [12] developed anisotropic and adjoint-based mesh adaptation for the DG discretization and cut-cell meshing. However, all these studies need output-based adjoint solutions as the weight in order to represent the sensitivity of the output functionals to the local primal residuals.
In shape optimization designs, the adjoint problems can be derived by two ways, i.e., the continuous approach and the discrete approach. Both of these approaches have made some successes. However, up to now, most of the adjoint-based adaptive algorithms have been developed only from the discrete adjoint approach [3][4][5][6][7][8][9][10][11][12][13][14][15] . In the discrete adjoint approach, discretization is made for the primal equations firstly, and then the adjoint equations are obtained by linearizing the discrete primal problems. In this procedure, a full Jacobian matrix is needed. Therefore, if we can get it from the primal flow solver, the implementation of the adjoint equation will be much easier. However, if not, the calculation of the full Jacobian matrix will be very complicated. Moreover, in the continuous adjoint approach, the adjoint partial differential equations (PDEs) are derived from the primal PDEs directly, and then these equations are discretized by the numerical schemes which may be different from the primal problem and may be independent of the code structure of the primal problem. More importantly, we can enforce the boundary conditions more rigorously in the continuous adjoint approach. Fidkowski and Roe [16] proposed an entropy adjoint approach, and used the continuous adjoint approach to derive the error indicators. Li et al. [17] developed an adaptive mesh refinement algorithm, and showed the promising results based on the continuous adjoint approach in the spectral difference framework. In this paper, we will provide an approach to construct the adaptive error indicator based on the continuous adjoint approach in the DG method, and show its efficiency for the error estimation.
The paper is organized as follows. In Section 2, we introduce the primal equations and the DG discretization used in this work. In Section 3, we introduce the derivation of the continuous adjoint equations and construct the adjoint-based error indicator. In Section 4, we use the continuous adjoint-based error to estimate the two-dimensional (2D) Euler equations. In Section 5, we describe the numerical implementation of our work. In Section 6, we discuss the numerical experiments to highlight the superiority of the proposed algorithm.

Governing equations
Let Ω be the computational domain in R 2 . Then, the 2D steady-state Euler equations can be written as follows: where u = (ρ, ρu, ρv, ρE) T is the solution vector consisting of the density ρ, V = (u, v) T is the velocity, and E is the total energy per unit mass. The flux functional F (u) is defined by where p is the pressure. The system is closed by where γ is the heat capacity ratio, which is 1.4 for the diatomic gas.

Standard DG discretization
Denote T h = {k} the subdivision of the domain Ω. Then, the standard DG discretization of Eq. (1) can be given as follows: find u h ∈ V h,p such that where In the above equations, d is the dimension of u, P p is the polynomial with the degree p, and (·) + and (·) − are the interior trace and the outer trace on ∂k, respectively.
, n) are the numerical flux functionals depending on the interior and ∂Ω meeting the physical boundary conditions, respectively.

Continuous adjoint equations
Consider the following primal problem arising from Eq. (1): where B is a boundary operator, Γ = ∂Ω, u ∈ V, and V is the functional space of the analytical solution. Then, the output functional which we are interested can be defined by where i Ω (·) and i Γ (·) are the operators which may be nonlinear. Let u h be the numerical approximation of u in the space V h,p ⊂ V. Then, we denote the error between u and u h by u, i.e., u = u h + u.
For the analytical solution u ∈ V and ψ ∈ V, the residual Therefore, we can add the residual to the output functional I(u) without changing its value, i.e., In this way, Therefore, the error of I(u) can be given by 1422 Huiqiang YUE, Tiegang LIU, and V. SHAYDUROV Then, we linearize the right-hand side of Eq. (10) to get the first-order approximations as follows: where the prime on (·) means the Fréchet derivative linearized at the term in the square brackets. Thus, we have In order to make I(u) − I(u h ) be independent of u, in other words, to find a suitable ψ such that the error u has no influence on I(u) − I(u h ), we consider the following problem: Therefore, ψ should satisfy Substituting Integrating the left-hand side of the above equation by parts, we get where the superscript * denotes the adjoint operation. Considering the right-hand side of Eq. (14), we can rewrite the continuous adjoint equations of Eq. (5) as follows:

Error estimation
In engineering calculation, we hope to get the best possible approximation of the output functionals with the smallest computing resources. An error estimate based on the computed values is required to evaluate how accurate these output functionals are approximated. From prior experience, it is very hard to estimate the error. Using the adjoint problem, we can obtain the general result.
Theorem 3.1 Let u be the analytical solution and u h be the numerical solution of the primal equation (5). Let ψ denote the analytical solution of the adjoint equation (16). Then, the error between I(u) and I(u h ) can be estimated by where Proof From this theorem, we can conclude that the primal residual can be related to the output functional error via the adjoint variable as a weight. Then, we take the absolute value of η k to construct the local error indicator for adaptive algorithms, namely, (20)

Some basic notations and output functionals
In this section, firstly, we define the flux Jacobian matrix as follows: Secondly, we define the boundary flux as follows: where n = (n 1 , n 2 ) is the normal vector of the boundary, and u n is the normal velocity, i.e., u n = un 1 + vn 2 . The pressure induced force coefficients, such as the drag coefficient, are the most important output functionals in inviscid compressible flows. Here, the drag coefficient is given by where Γ = ∂Ω, and Γ W is the solid wall of Γ. C ∞ = 1 2 ρ ∞ |V ∞ | 2 l, where l is the reference length. ϕ = (cos(θ), sin(θ)) T , where θ denotes the angle of attack. The subscript ∞ means the far-field quantity.

Adjoint equation and boundary conditions
In this section, based on the drag coefficient, we introduce how to derive the continuous adjoint Euler equations and its boundary conditions. The approaches for other output functionals are similar. From Eq. (21), we note that the drag coefficient does not depend on the domain terms, i.e., i Ω (u) = 0 in Eq. (7). Therefore, the adjoint equation (16) can be given by where ψ = (ψ 1 , ψ 2 , ψ 3 , ψ 4 ) is the adjoint solution vector.
In inviscid compressible flows, on the solid wall Γ W , the boundary condition imposes u n = 0. Therefore, n · F(u) = p(0, n 1 , n 2 , 0) T on Γ W . Recalling the general functional of the boundary condition (17), the adjoint boundary condition on Γ W under this output functional can be given by i.e., On a far-field boundary ∂Ω ∞ , recall the definition of the drag coefficient (21). We notice that the value of this output functional does not depend on the far-field. Therefore, simply, we set Finally, we obtain the adjoint problem as follows: 5 Numerical implementation

Solving primal equations
Firstly, we add a pseudo time term to Eq. (1), i.e., we rewrite the steady-state Euler equations as follows: The spatial discretization by the DG method leads to a system of ordinary differential equations as follows: where M is the mass matrix, and R is the right-hand residual. In this work, we are interested in stationary problems. Therefore, the backward Euler scheme is used in the temporal discretization. After temporal discretization and linearizing Eq. (28) in time, we can get a system of linear equations as follows: where Δt is the time increment, and Then, the linear system can be solved by the generalized minimal residual (GMRES) method with ILU0 preconditioning.

Solving adjoint equations
In Eq. (19), the error indicator involves the analytical adjoint solution ψ which is unknown. Therefore, we need to make an approximation to it. Let ψ h denote the approximation. We note that ψ h must not be approximated by the same finite element space V h,p . Otherwise, the error indicator will be identically equal to zero due to the Galerkin orthogonality. There are several approaches to compute the numerical approximation ψ h [14][15] . In this work, the adjoint equations are approximated on the same mesh employed for the primal equations, but the order of the finite element space is increased from p to p + 1.
The numerical implementation of the adjoint equations is a little different from the primal equations. From Eq. (26), we notice that (A T 1 , A T 2 ) is independent of the adjoint variable ψ. Therefore, the adjoint equations are linear. Therefore, we can solve this problem by iteration or the direct linear system solver directly after spatial discretization.
We should remark that there are two main differences between the DG discretization of the primal problem and the DG discretization of the adjoin problem. The first one is the difference between the numerical fluxes. Here, we write Eq. (1) in the following quasi-linear form: i.e., (A 1 , A 2 ) · ∇u = 0 in Ω. We note that Eq. (22) and the above equation have the same eigenvalues. Therefore, the local Lax-Friedrichs numerical flux for the adjoint equation can be given by where α is the largest eigenvalue which can be calculated from the primal problem directly. The second difference lies in the boundary conditions. As we have discussed in Section 4 on the far-field boundary ∂Ω ∞ , we simply set ψ − h = 0. On the solid wall boundary, there are a set of boundary values of ψ h satisfying the boundary conditions in Eq. (26). Here, we enforce the solid wall boundary condition as follows: which guarantees the average value satisfying the adjoint solid wall boundary condition, i.e., ψ 2 n 1 + ψ 3 n 2 = 1 C∞ n · ϕ.

Adjoint-based adaptive mesh refinement algorithm
Utilize the local error indicator (20) by an adaptive mesh refinement algorithm shown in Fig. 1. In each refinement cycle, the cells with a large local error indicator η ind k will be refined, and the cells with small η ind k will be coarsened. Following the process, we will get a mesh on which the local error indicator η ind k is equidistribution. For more details of the mesh adaptation strategies, we refer readers to Ref. [8].

Numerical results
We show some numerical tests to discuss the advantages of our adaptive mesh refinement algorithm, comparing with the traditional algorithms based on the feature-based or residualbased error indicators. Here, firstly, let us introduce the feature-based error indicator η feature k = lg(1 + |∇ρ k |) and the residual-based error indicator [18] η residual where ∇ρ k is the gradient of the density, and In the following numerical tests, the primal problem will be approximated in V h,1 , and the adjoint problem will be approximated in V h,2 .

NACA-0012 airfoil
In the first numerical experiment, we consider an NACA-0012 airfoil in an inviscid compressible flow with M a ∞ = 0.4 and θ = 0. The output functional considered here is the drag coefficient. Therefore, the theoretical value of the output functional is zero. Figure 2 shows the initial mesh of this test. The meshes obtained by use of both the adjoint error indicator and the feature-based error indicator are shown in Fig. 3. From these figures, we can see that most of the refinement is enforced near the leading and trailing edges of the airfoil. Furthermore, we find that the mesh obtained by our adjoint-based error indicator is also refined around the surfaces of the airfoil. Figure 4 shows the output functional error convergence under different refinement approaches. From these results, we can see that the error of the output functional estimated on the mesh obtained by use of the adjoint-based error indicator is the smallest for a given number of cells. In addition, we can observe that, at a given error, the number of the mesh cells needed by the adjoint-based adaptive approach is an order of magnitude less than the uniform refinement approach. Figure 5 shows the distribution of the adjoint-based error indicator during the refinement cycle. We find that the errors are larger around the surfaces of the airfoil on the initial mesh.   The cells with larger errors are refined, which explains why the cells near the airfoil surfaces are refined. Also, we notice that with the repetition of the mesh refinement, the errors in all cells are decreasing. Consequently, using our adaptive algorithm, we can get a mesh, on which the error is equally distributed.

Smooth bump
In the second numerical experiment, we consider the inviscid flow over a smooth bump. The computational domain is given by Ω = {(x, y) ∈ [−1.5, 1.5] × [0, 0.8], y > 0.062 5 e −25x 2 }. At the inlet, M a ∞ = 0.5, and θ = 0. There should be no entropy production because the flow is smooth enough. Therefore, we construct the following output functional: where p is the pressure, ρ is the density, and the subscript ∞ means the freestream value. At the left-and right-sides, we apply the subsonic inflow boundary condition and the subsonic outflow boundary condition, respectively. At the top and the bottom, the slip wall conditions are applied. The theoretical value of this output functional is zero. Figure 6 shows the initial mesh for this problem, and Fig. 7 shows the pressure contour of the primal Euler equations. We can see that the solution is axisymmetric. From Fig. 8, we can see that most of the refinements are near the bump. Both of the two adaptive approaches construct almost the same meshes, except that the refinement is denser around the bump in the adjoint-based approach.
In Fig. 9, we show the convergence of the output functional error. Some observation can be made. Firstly, both the adjoint-based and residual-based adaptive approaches show better convergence than the uniform refinement. Indeed, after several refinements, the error computed on the subsequent adaptive mesh is at least two orders of magnitude smaller comparing with the error computed on the uniform refinement mesh. Secondly, the adjoint-based refinement gives a better accuracy than the residual-based refinement for a given number of cells.

Conclusions
We have introduced an adjoint-based adaptive mesh refinement algorithm for inviscid compressible flows. Firstly, we show the derivation of a continuous adjoint problem for the Euler equations. Secondly, we introduce the error estimation for the given output functionals. Thirdly, we discuss the numerical implementation of the DG method for an adjoint problem. Then, we develop an adaptive mesh refinement algorithm based on the continuous adjoint approach. The numerical tests show that the meshes constructed by this adaptive algorithm are much more economical for computing the values of the output functionals than the meshes constructed by the traditional error indicators. We conclude that the adaptive approach presented in this paper is very promising in improving the accuracy of the output functionals within given computational resources.